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1.  Introduction 


Under  Office  of  Naval  Research  sponsorship,  the  Marine  Physical  Laboratory  has 
designed,  built  and  deployed  a  set  of  12  acoustic  sensors  which  are  neutrally  buoyant  and  drift 
freely  when  deployed  in  the  ocean.  The  sensors  are  called  Swallow  floats  after  J.  C.  Swallow 
who  used  freely-drifted  floats  to  measure  deep  ocean  currents.*  The  MPL  Swallow  floats  also 
drift  freely,  but  are  used  to  measure  acoustic  energy  in  the  very  low  frequency  (VLF)  band, 
which  is  generally  defined  as  1  to  20  Hz. 

The  Swallow  floats  are  designed  to  operate  without  tethers  so  that  their  measurements 
are  not  contaminated  by  tether  strumming  noise  and  so  that  they  drift  with  the  ocean  current 
and  are  not  subject  to  flow  noise.  Each  float  possesses  an  8  kHz  acoustic  transducer  for  measur¬ 
ing  interfloat  and  float-to-surface  travel  times.  The  travel  times  may  be  used  to  estimate  time- 
varying.  relative  float  positions.  Several  floats  deployed  in  the  same  general  area  form  a  freely- 
drifting,  randomly-spaced  array  whose  acoustic  measurements  may  be  combined,  or  bearnformed. 

The  MPL  Swallow  floats  have  been  deployed  several  times  annually  since  1982  as  part 
of  the  system  development  process.  In  terms  of  analyzing  and  interpreting  Swallow  float  meas¬ 
urements,  single  float  measurements  of  ambient  ocean  noise  have  been  shown  to  be  consistent 
with  other,  similar  measurements  found  in  the  literature.*  This  thesis  describes  how  float  posi¬ 
tions  are  estimated  from  interfloat  travel  time  measurements,  and  demonstrates  the  correct  way 
to  beamform  signals  from  the  directional  geophones.  The  localization  and  beamforming  pro¬ 
cedures  are  illustrated  using  measurements  from  a  September  1986  Swallow  float  experiment 
conducted  approximately  50  miles  west  of  San  Diego,  California. 


2.  Swallow  Float  System  Description 

A  number  of  articles  have  been  published  to  date  providing  details  of  the  Swallow  float 
hardware  and  operation. Reference  [6|  is  probably  the  most  complete  system  description. 
The  intent  here  is  to  provide  sufficient  detail  to  enable  the  reader  to  understand  the  results 
which  follow. 

As  illustrated  in  Figure  2.1,  the  Swallow  float  system  consists  of  a  17  inch  (0.432  meter) 
diameter,  9/16  inch  (11.3  mm)  thick  glass  sphere  which  contains 

•  3  orthogonally-oriented  velocity  sensors,  known  as  geophones; 

•  a  control  [irocessor  unit  and  electronic  circuitry  to  condition,  sample  and 
record  the  geophone  signal,  and  to  control  other  float  activities  such  as  ballast 
release; 

•  a  compass  which  measures  the  direction  of  the  y  axis  geophone  relative  to 
magnetic  north; 

•  a  solid  state  memory  buffer  and  a  digital  cartridge  tape  recorder; 

•  a  battery  pack  power  supply. 

External  to  the  sphere  are 

•  an  8  kHz  acoustic  transducer  for  transmitting  and  receiving  ranging  signals; 

•  an  optical  flasher  and  a  radio  beacon  used  to  locate  the  float  when  it  is  on  the 
surface; 
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Figure  2.1  Swallow  float  hardware 
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•  expendable  ballast  and  a  ballast  release  mechanism; 

•  a  plastic  hardhat  which  protects  the  glass  sphere. 

In  operation  the  floats  are  lowered  into  the  water  from  the  deck  of  a  ship  after  which 
they  sink  to  a  preset  depth  and  operate  autonomously.  They  are  ballasted  to  neutral  buoyancy 
at  a  particular  depth,  or  in  some  cases,  to  sink  to  the  ocean  bottom.  Recent  experiments  have 
utilized  3  floats  positioned  on  the  ocean  bottom  to  pi'ovide  a  stationary  coordinate  system 
within  which  the  remaining  floats  can  be  localized.  Upon  command  from  the  surface,  or  failing 
to  receive  that  command,  according  to  a  preset  timer,  the  floats  release  a  portion  of  their  bal¬ 
last  and  ascend  to  the  surface. 

While  deployed,  each  float  samples  and  internally  records  signals  from  the  3  geophones, 
the  compass,  and  the  8  kHz  range  pulse  detection  circuit.  The  floats  divide  time  into  45  second 
records  consisting  of  44  seconds  during  which  data  are  sampled  and  stored  in  the  buffer  and  1 
second  during  which  the  buffer  contents  are  written  to  tape.  No  geophone  data  are  sampled 
during  the  1  second  buffer  dump.  Appendix  1  of  Reference  [13]  contains  a  detailed  description  of 
the  conditioning  applied  to  the  geophone  signal.  Briefly,  the  geophones  produce  a  signal  propor¬ 
tional  to  their  velocity,  which  is  then  amplified,  low-pass  filtered,  and  sampled  at  50  Hz.  Float 
data  storage  capacity  is  limited  to  17  Mbytes,  which  corresponds  to  about  25  hours  of  opera¬ 
tion. 

Once  a  float  has  surfaced  and  been  retrieved,  its  data  tape  is  extracted  and  read.  The 
8  kHz  acoustic  pulse  arrival  times  provide  interfloat  and  float-to-surface  travel  time  measure¬ 
ments  which  may  be  used  to  determine  the  float  positions.  Geophone  signals  provide  float  velo¬ 
city  measurements  from  which  acoustic  pressure  may  be  estimated.  Float  compass  heading 
measurements  arc  made  because  the  geophones  are  directional  sensors  and  we  must  know  the 
direction  in  which  they  were  pointed  when  the  measurements  were  made. 


3.  September  1986  Experiraent 

Data  presented  in  this  thesis  were  acquired  between  16  and  19  September  1986, 
approximately  50  miles  west  of  San  Diego,  California  and  approximately  20  miles  southeast  of 
San  Clemente  Island  (32  °  35' N,  118' 10'  W).  Figure  3.1  shows  the  location  of  the  deployment 
site.  A  detailed  look  at  individual  float  data  from  this  experiment  is  contained  in  Reference  [9]. 
This  Chapter  contains  a  summary  of  (9)  and  some  additional  information  necessary  to  localize 
the  floats  and  beamform  their  outputs. 


3.1.  Experiment  Summary 

The  experiment  began  on  Tuesday,  16  September  1986,  when  6  MPL  personnel  with  12 
Swallow  floats  sailed  from  the  University  of  California’s  Nimitz  Marine  Facility  in  San  Diego 
aboard  the  110'  chartered  vessel  Scorpiua.  The  experiment  site  was  the  northwest  corner  of  San 
Clemente  Basin,  where  the  water  depth  is  about  1000  fathoms  (about  2000  meters).  Three  of 
the  floats  (0,  1,  and  2)  were  ballasted  to  descend  to  the  ocean  bottom.  The  remainder  were  bal¬ 
lasted  for  neutral  buoyancy  at  1000  meters  depth.  The  experiment  was  scheduled  to  begin  on 
the  morning  of  17  September  and  to  last  about  24  hours.  Inside  each  buoy,  the  backup  ballast 
release  timers  had  been  set  to  release  on  19  September. 

Figure  3.2  shows  the  bottom  elevation  contours  in  the  vicinity  of  the  experiment  site. 
San  Clemente  Basin  is  bounded  to  the  east  by  a  series  of  underwater  cliffs  which  lead  up  to  very 
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shallow  water,  less  that  100  meters  deep  in  some  places.  The  Basin  is  bounded  to  the  north  by 
upward  sloping  ground  which  rises  to  San  Clemente  Island.  To  the  west  and  south,  ihe  bottom 
slopes  upward  more  gradually. 

The  Scorpius  arrived  on  station  early  on  the  17^.  The  floats  were  checked  out  and 
found  to  be  operating  satisfactorily.  They  were  initialized  synchronously  at  0515,  and  by  1000 
all  12  floats  had  been  deployed.  Floats  6,  7,  8  and  11  encountered  problems  and  had  to  be 
dropped  from  the  experiment,  Figure  3.3  (taken  from  Figure  1.3  of  [9])  shows  where  the  floats 
were  put  in  the  water  and  later  recovered. 

Scorpius  remained  in  tiic  general  vicinity  of  the  experiment  site  throughout  the  day  on 
the  17*^,  and  in  the  evening  sailed  to  .San  Clemente  Island,  where  she  was  anchored  for  the 
night.  The  ship  got  underway  again  at  0100  on  the  18^  and  arrived  at  the  experiment  site  at 
about  0630.  The  day  was  spent  recalling  and  recovering  floats  one  at  a  time.  By  1600  the  last 
float  was  aboard  and  the  ship  sailed  for  San  Diego. 


3.2.  Acoustic  Range  Data 

Acoustic  range  pulses  transmitted  by  a  float  are  reflected  by  the  surface  and  received 
by  that  same  float.  Figure  3.4,  taken  from  Figure  111.12  of  (9),  shows  float  lO’s  surface  echoes. 
The  vertical  axis  has  been  scaled  from  round  trip  travel  time  to  depth  using  one-half  the  sound 
speed,  approximately  750  metcrs/second.  The  horizontal  axis  is  elapsed  deployment  time.  The 
echoes  from  the  surface  descend  gradually  from  the  left  and  level  off  at  about  1000  meters 

depth.  The  echoes  from  the  ocean  bottom  rise  from  the  left  and  obscure  the  surface  bounce 

pulse  arrival  between  records  1000  and  1(50.  This  establishes  record  1150  as  the  first  record 
where  we  can  begin  to  localize  the  floats.  Table  3.1  shows  the  depths  to  which  the  floats 
deployed  based  upon  their  surface  echo  data. 

Figure  3.5,  taken  from  Figure  IV.lc  of  [9],  shows  float  O’s  record  of  pulses  received  from 
float  3  (freely-drifting  float  to  bottomed  float).  The  short  dashes  extending  across  the  page  are 
the  direct  path  pulse  (10  msec,  or  150  meters  long)  and  the  longer  dashes  are  the  surface  bounce 
arrivals.  Scattering  by  the  rough  sea  surface  is  thought  to  be  responsible  for  the  extended  dura¬ 
tion  of  the  surface  bounce  arrivals.  The  direct  and  surface  bounce  arrivals  come  together 
quickly  beginning  about  record  1900,  marking  float  3’s  ballast  release  and  ascent  to  the  surface 
and  also  indicating  the  end  of  the  period  during  which  we  may  localize  the  floats. 

Direct  path  acoustic  range  pulses  arc  clearly  visible  between  all  pairs  of  floats  except 

between  the  bottomed  floats.  The  upward  refracting  sound  speed  gradient  apparently  caused 

severe  attenuation  of  direct  path  transmission  between  the  bottomed  floats.  However,  the  sur¬ 
face  bounce  path  is  clearly  visible  for  these  floats.  Figure  3.6,  taken  from  Figure  FV.lb  of  [9], 
shows  the  arrival  of  pulses  transmitted  by  float  2,  reflected  by  the  surface,  and  received  by  float 
0. 

Distinct  pulse  arrival  times  are  seen  for  all  pairs  of  floats  and  the  surface  echoes 
between  records  1150  and  1870.  The  po.sitions  of  the  floats  during  this  interval  will  be 
estimated  in  Chapter  4. 


3.3.  AGO  Level  and  Float  Heading 

Reference  [9]  contains  plots  of  automatic  gain  control  (AGC)  level  and  float  compass 
heading  for  each  float  during  the  deployment.  One  such  plot,  for  float  3,  is  presented  in  Figure 
3.7.  The  AGC  indicates  the  amount  of  gain  applied  to  the  geophone  signals.  It  is  programmed 
to  increase  as  signal  decreases,  and  vice  versa,  in  order  to  optimize  the  dynamic  range  of  the 
analog  to  digital  converter  (ADC).  The  AGC  level  is  thus  inversely  proportional  to,  and  a 
rough  measure  of,  geophone  signal  level.  In  Figure  3.7,  the  AGC  level  is  zero  during  the  float’s 
descent  because  of  the  descent  motion,  and  higher  thereafter  once  the  float’s  depth  becomes 
stable.  The  drop  in  AGC  level  between  records  1400  and  1600  was  caused  by  high  level  signals 
received  from  a  pa.s.sing  ship. 
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Table  3.1:  Deployment  Depths 

Float  Denth.  meters 

0 

2010 

(bottom) 

1 

1820 

(bottom) 

2 

1900 

(bottom) 

3 

960 

4 

1040 

5 

1060 

9 

800 

10 

990 

Float  compass  heading  is  the  direction  of  the  y  axis  geophone,  measured  clockwise  from 
magnetic  north.  The  compass  heading  trace  in  Figure  3.7  shows  that  float  3  was  on  the  deck  of 
the  ship  until  about  record  200,  that  it  rotated  as  it  descended  between  records  200  and  750, 
and  that  its  depth  (and  to  a  large  degree,  its  heading)  stabilized  after  about  record  750. 


3.4.  Geophone  Data 

As  a  first  step  in  evaluating  Swallow  float  geophone  data,  the  rms  power  in  5  second 
segments  is  computed.  The  resulting  power  series,  which  are  shown  in  [9],  show  gross  signal 
characteristics.  One  such  plot,  for  records  1750  through  1980  from  float  5,  is  shown  in  Figure 
3.8.  The  power  series  showed  that  three  of  the  floats  (0,  2  and  4)  had  connector  problems  which 
grounded  one  geophone  channel.  This  was  not  a  problem  in  the  case  of  floats  0  and  2,  since 
they  were  tethered  to  the  bottom  and  their  geophone  signals  are  almost  entirely  clipped  due  to 
motion  from  tether  strumming.  For  float  4,  however,  it  meant  that  the  float  could  not  be  used 
in  the  beamforming. 

A  problem  discussed  in  [9]  is  the  presence  of  strong  oscillatory  signals  at  the  start  of 
each  record  which  decay  exponentially.  These  signals  are  attributed  to  mechanical  oscillation 
of  the  float,  apparently  excited  by  the  internal  tape  recorder  when  it  cycles  on  and  off  between 
records.  They  can  cause  velocity  data  to  be  clipped  for  up  to  8  seconds  at  the  beginning  of  the 
record,  obscuring  the  ambient  ocean  noise  signal  and  reducing  the  dynamic  range  with  which 
the  ambient  ocean  noise  is  measured  during  the  remainder  of  the  record. 

Since  the  September  1986  deployment,  the  Swallow  floats  have  been  modified  to  elim¬ 
inate  or  greatly  attenuate  mechanical  vibrations.  It  was  also  desirable  to  remove  these  signals 
from  the  1986  data  set  to  permit  its  further  use.  Reference  [lO]  derives  a  method  of  removing  or 
significantly  reducing  resonant  oscillation  present  in  data,  which  may  be  applied  if  certain  con¬ 
ditions  are  met.  The  method  is  shown  to  be  effective  on  the  1986  deployment  data  set  and  was 
applied  to  data  from  floats  3,  9,  and  10  before  continuing  with  the  analysis. 

A  related  problem,  which  is  addressed  separately  in  Reference  (ll],  is  resonant  float 
oscillation  at  frequencies  below  0.6  Hz.  Oscillation  at  these  frequencies  appears  to  consist  of  an 
irregular  rocking  or  wobble  which  is  strongly  excited  as  the  floats  descend,  weakly  but  con¬ 
sistently  excited  by  the  tape  recorder  at  the  start  of  each  record,  and  occasionally  excited  by 
the  receipt  of  strong  acoustic  pulses.  We  are  not  directly  concerned  with  this  problem  here 
because  it  is  below  the  band  of  immediate  interest. 

The  focus  of  the  beamforming  in  Chapter  5  will  be  on  records  1760  through  1870, 
which  is  the  last  110  records,  or  82.5  minutes,  of  the  period  during  which  the  float  positions  will 
be  estimated  in  Chapter  4.  This  portion  of  the  data  record  was  chosen  with  the  hope  of  seeing 
propulsion  tonals  from  Seorpius  as  she  steamed  down  from  San  Clemente  Island  early  on  the 
morning  of  18  September.  Since  both  engines  were  running  at  1750  rpm  and  the  reduction  gear 
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ratio  is  known  to  be  3.8:1,  the  shaft  rate  may  be  expected  to  be 

'If  .  ' 

/sVaft  ra  <■  =  (^^^0  r«v./min.)  •  - - - -  •  =7.67  Hz  (3.1) 

■'  shaft  race  \  &  /  '  [  3.8  engine  rev./shafl  rev.  j  60  sec  J  ^  ' 

The  propeller  has  4  blades  to  that  the  blade  pass  rate  is  30.7  Hz.  which  is  out  of  band. 

Appendix  A  contains  plots  of  some  of  the  time  series  from  the  selected  period  and 
Appendix  B  contains  the  corresponding  spectra.  Data  from  floats  3,  9  and  10  have  been  pro¬ 
cessed  as  described  in  Reference  [lO]  to  remove  resonant  mechanical  oscillation  signals  excited 
by  the  tape  recorder  turning  on  and  off.  A  2048  point  (40.96  second)  FFT  was  used  in  Appendix 
B,  and  data  were  windowed  with  a  Kaiser-Bessel  window  (a  =  2.5).  Only  positive  frequencies 
are  shown.  Reference  [2]  discussed  calibrating  the  spectra. 

Several  spectral  lines  are  visible  in  all  the  records.  A  tone  is  seen  at  6.7  Hz  in  some 
records,  which  may  be  Scorpius’s  shaft  rate,  but  there  are  other  spectral  lines  at  much  higher 
levels.  A  set  of  lines  at  9,  13.5,  18  and  22.5  Hz  are  harmonically  related,  although  the  funda¬ 
mental  frequency,  4.5  Hz,  is  absent.  Two  other  lines,  at  10.3  and  12.2  Hz,  are  also  prominent. 

Figure  3.9  shows  how  the  frequencies  of  these  lines  vary  with  time.  The  plots  were 
made  by  searching  a  small  part  of  the  spectrum  near  each  frequency  for  the  exact  bin  at  which 
the  peak  occurs.  Note  that  the  frequency  of  the  lines  at  approximately  9  and  13.5  Hz  slowly 
decreases  during  the  HO  record  interval.  All  the  floats  record  rather  abrupt  shifts  in  these  fre¬ 
quencies  near  records  1850  and  1865,  and  in  the  9  Hz  line  only,  near  record  1808.  These  abrupt 
frequency  shifts  also  occur  in  the  lines  at  approximately  18  Hz  and  22.5  Hz,  indicating  that  the 
9,  13.5,  18  and  22.5  Hz  lines  are  all  harmonically  related.  The  10.3  and  12.2  Hz  lines  are  stable 
throughout  the  interval,  indicating  that  they  are  not  harmonically  related  to  the  other  lines. 

The  source  of  all  of  the  lines  is  assumed  to  be  a  single  surface  ship.  It  will  become 
apparent  in  Chapter  5  that  all  of  the  lines  arrive  from  approximately  the  same  bearing.  The  9 
Hz  and  harmonics  line  set  is  thought  to  be  propulsion  related  because  the  fundamental  fre¬ 
quency  (4.5  Hz)  is  normal  for  a  merchant  ship  blade  pass  frequency  and  is  usually  accompanied 
by  strong  harmonics.*®  The  reason  for  the  absence  of  the  fundamental  is  hypothesized  in 
Chapter  6.  The  10.3  and  12.2  Hz  lines  are  likely  to  be  auxiliary  machinery  tones. 

Figure  3  10  shows  the  level  of  the  9  Hz  line,  which  rises  by  about  10  dB  during  the 
intervil.  The  levels  of  the  other  lines  are  not  shown  here,  but  they  decrease  or  remain  approxi¬ 
mately  constant  during  the  interval.  Thus  nothing  about  the  range  to  the  source  can  be  con¬ 
cluded.  However,  the  transmission  losses  are  apparently  different  for  different  frequencies 
emanating  from  the  same  source. 

Figure  3.11  shows  the  vertical  angle  of  arrival  for  the  9,  10.3,  12.2,  and  13.5  Hz  signals 
during  the  interval.  The  arrival  angle  was  obtained  by  comparing  the  horizontal  and  vertical 
spectral  levels.  For  a  particular  frequency  and  during  any  one  record,  the  arrival  angle  is 
different  at  each  float.  Also,  for  a  particular  float  and  during  any  one  record,  the  arrival  angle 
is  different  at  each  frequency.  This  indicates  that  the  sound  field  is  far  from  homogeneous,  and 
that  sound  is  not  propagating  across  the  Swallow  float  array  as  plane  waves.  The  beamforming 
results  of  Chapter  5  will  be  more  understandable  in  light  of  Figure  3.11. 


3.5.  Sound  Speed  Profile 

The  sound  speed  profile  for  the  upper  500  meters  was  obtained  from  temperature  meas¬ 
urements  made  by  the  California  Cooperative  Oceanic  Fisheries  Investigations  (CalCOFI)  Cruise 
8609  on  19  September  1986  at  (32  ”  30.9'  N,  1 18  ”13.3'  W),*^  which  is  within  a  few  miles  of  the 
deployment  site.  Historical  data  archived  by  the  National  Oceanographic  Data  Center  were 
used  for  salinity  in  the  upper  500  meters  and  temperature  and  salinity  in  the  lower  1500 
meters.*®  Using  the  equation  given  by  Mackenzie,*®  the  composite  sound  speed  profile  shown  in 
Figure  3.12  was  calculated. 
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Figure  3.9a  Exact  frequency  of  spectral  lines  near  9  Hz  (upper)  and 
13.5  Hz  (lower  panel) 
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Figure  3.9b  Exact  frc(|ucncy  ol  spectral  lines  near  10.3  Hz  (upper)  and 
12.2  Hz  (lower  panel) 
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Figure  3.10  Spectral  level  of  the  line  near  9  Hz 


4.  Swallow  Float  Localization 


In  this  Chapter,  the  time-varying  locations  of  the  Swallow  floats  are  estimated  for 
records  1150  to  1870  of  the  September  1986  experiment.  The  rms  position  error  is  also 
estimated.  The  goal  of  the  localization  is  to  obtain  a  position  estimate  with  an  rms  error  of  less 
than  A/10  at  the  highest  frequency  of  interest,  since  theory  predicts  that  beamforming  results 
will  be  degraded  by  less  than  1  dB  if  this  condition  is  met.*®  Since  the  highest  frequency  of 
interest  is  20  Hz,  the  condition  to  be  satisfied  is 


rms  position  error  < 


X  _  c//  _  1500  m/sec 
10  “  10  ”  10(20  Hz) 


7.5  meters 


(4.1) 


4.1.  Theoretical  Development  and  Simulation  Results 

Aspects  of  the  localization  problem  are  addressed  in  an  MPL  Technical  Memorandum^* 
and  in  a  Journal  of  Oceanic  Engineering  (JOE)  submission,*^  the  latter  of  which  is  included  here 
as  Appendix  C.  The  earlier  results  are  summarized  below  and  extended  in  subsequent  subsec¬ 
tions. 
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Figure  3.11b  Magnitude  of  the  vertical  arrival  angle  at  10.3  Hz  (upper) 
and  12.2  Hz  (lower  panel) 
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4.1.1.  Summary  of  Previous  Results 

Reference  [21]  describes  how  interfloat  travel  times  and  travel  time  error  variance  are 
obtained  from  the  acoustic  range  data.  Briefly,  the  pulse  arrival  times  are  detected,  and  the 
pulse  transmit  time  subtracted  from  them,  to  obtain  measurements  of  the  travel  time  between 
each  float  pair  and  between  each  float  and  the  surface.  Only  one  float  transmits  during  each 
record,  so  a  particular  travel  time  measurement  is  made  at  N  record  intervals,  where  N  is  the 
number  of  floats.  The  travel  time  measurements  must  be  interpolated  to  obtain  a  measurement 
for  each  record.  Reciprocal  path  measurements  are  then  averaged  to  obtain  travel  time  esti¬ 
mates.  The  difference  between  the  reciprocal  path  measurements  is  used  to  estimate  the  travel 
time  error,  as  shown  in  [21].  A  diagram  showing  the  processing  steps  may  be  found  in  Figure  4 
of  Appendix  C. 

Two  methods  are  derived  in  Appendix  C  for  estimating  float  positions  from  travel  time 
and  travel  time  error  estimates: 

•  A  generalized  least  squares  filter  which  estimates  current  position  using  the  current 
set  of  travel  time  measurements.  The  filter  uses  only  that  set  of  measurements 
and  an  estimate  of  the  rms  measurement  error. 

•  A  Kalman  filter  which  estimates  current  position  by  weighting  the  current  set  of 
measurements  against  the  previous  position  estimate  propagated  forward  in  time. 
Position  estimates  are  propagated  forward  in  time  using  equations  of  motion  for 
the  floats.  The  filter  must  also  be  supplied  with  estimates  of  the  rms  measurement 
error  and  float  accelerations. 

Simulations  are  employed  in  Appendix  C  to  test  the  filters  under  a  variety  of  deploy¬ 
ment  geometries  and  noise  conditions.  Four  simulation  deployment  geometries  are  shown  in  Fig¬ 
ure  6  of  Appendix  C.  Filter  performance  is  evaluated  for  four  combinations  of  process  and 
measurement  noise.  Process  noise  parameterizes  how  erratically  the  floats  move  between  travel 
time  measurements,  and  measurement  noise  is  the  error  in  the  travel  time  measurements.  For 
the  simulations,  float  accelerations  were  modeled  as  normally  distributed  white  noise,  and  pro¬ 
cess  noise  is  actually  the  variance  of  the  acceleration  used  in  the  simulation. 

Cases  1,2,6  and  7  test  the  most  frequently  used  Swallow  float  deployment  geometry 
in  which  3  floats  are  deployed  to  the  ocean  bottom  and  the  remainder  in  a  planar  array  some¬ 
where  in  the  water  column.  Case  1  uses  low  process  and  measurement  noise  in  order  to  estab¬ 
lish  the  baseline  for  evaluating  how  the  filters  respond  to  increased  noise  or  other  deployment 
geometries.  Cases  2,  6,  and  7  test  the  filters’  performance  when  the  process  and  measurement 
noise  are  increased.  The  rms  measurement  noise  is  raised  by  a  factor  of  2  in  Cases  2  and  7,  and 
the  rms  process  noise  is  raised  by  a  factor  of  about  100  in  Cases  6  and  7. 

Both  the  process  and  measurement  noise  are  low  in  Cases  3,  4  and  5.  Case  3  tests  the 
filters’  performance  when  the  mid-water  column  floats  drift  out  of  the  net  formed  by  the  bot¬ 
tomed  floats.  The  bottomed  floats  are  the  same  as  in  Case  1,  but  the  planar  array  is  now  out¬ 
side  the  triangle  formed  by  the  bottomed  floats.  Case  4  tests  the  effect  of  one  of  the  bottomed 
floats  failing.  The  remaining  bottomed  floats  and  planar  array  are  the  same  as  in  Case  1.  Case 
5  tests  the  effect  of  deploying  the  mid-water  column  floats  in  a  vertical  line  array  rather  than  a 
planar  array.  The  three  bottomed  floats  were  put  at  about  4000  meters  depth  and  5  other 
floats  spaced  between  700  and  2200  meters  depth.  The  simulation  cases  are  summarized  in 
Table  4.1. 

The  coordinate  system  for  the  localization  is  shown  in  Figure  4.1.  The  origin  is  fixed  at 
the  surface  directly  above  one  of  the  bottomed  floats.  This  fixes  the  location  of  the  Z  axis.  The 
horizontal  axes  are  then  rotated  until  a  second  bottomed  float  lies  in  the  XZ  plane. 

The  results  of  the  simulations  are  shown  in  Figures  7  and  8  of  Appendix  C.  Figure  7 
compares  the  results  of  Cases  1,  2,  6  and  7  and  shows  the  effects  of  increased  process  and  meas¬ 
urement  noise.  The  horizontal  axis  is  process  noise  and  the  vertical  axis  is  rms  position  error. 
The  two  lower  lines  show  the  results  when  the  measurement  noise  is  low,  and  the  upper  two 
when  the  measurement  noise  is  high.  The  dotted  lines,  which  depict  the  least  squares  filter 
results,  are  parallel  to  the  horizontal  axis  because  the  filter  is  insensitive  to  the  level  of  process 
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Table  4.1.  Float  Location  Simulations. 

Case 

Process  Noise 

Measurement  Noise 

Tests  the  Effect  of; 

1 

low 

low 

filter  performance  baseline 

2 

low 

high 

increased  measurement  noise 

7 

high 

low 

increased  process  noise 

6 

high 

high 

increased  process  and  measurement  noise 

3 

low 

low 

large  array  drift 

4 

low 

low 

losing  a  bottomed  float 

5 

low 

low 

vertical  line  array 

noise,  The  filter  processes  one  set  of  measurements  at  a  time  and  is  not  affected  by  what  hap¬ 
pens  to  the  floats  between  measurements. 

We  see  that  the  error  in  the  Kalman  filter  estimate  is  smaller  than  that  in  the  least 
squares  estimate  as  long  as  the  process  noise  is  low.  This  is  because  the  Kalman  filter  tracks 
the  floats’  motion  when  the  process  noise  is  low.  As  the  process  noise  is  increased,  the  Kalman 
filter  is  less  able  to  maintain  a  track  and  the  error  in  its  position  estimate  increases.  At 
sufficiently  high  process  noise,  the  two  filters  perform  comparably.  The  effect  of  increasing  the 
measurement  error  is  to  increase  the  error  in  the  output  of  both  filters. 

Figure  8  compares  the  results  of  Cases  1,  3,  4  and  5  which  show  the  effect  of  different 
float  deployment  geometries  on  the  filters’  performance.  The  Kalman  filter  performs  somewhat 
better  than  the  least  squares  filter  in  all  of  the  deployment  geometries  simulated.  Since  these 
simulations  were  conducted  with  low  process  noise,  the  performance  differential  is  probably  due 
to  the  tracking  capability  of  the  Kalman  filter. 

There  remain  three  issues  which  must  be  addressed  prior  to  applying  the  Kalman  and 
least  squares  filters  to  experimental  data.  They  are: 

•  obtaining  a  good  initial  estimate  of  float  positions  and  velocities, 

•  determining  the  relationship  between  actual  position  error  and  the  error  esti¬ 
mates  produced  by  the  filters; 

•  determining  how  error  in  the  sound  speed  estimate  impacts  the  performance  of 
both  filters; 

These  issues  are  discussed  below. 


4.1.2.  Obtaining  Good  Initial  Position  and  Velocity  Estimates 

Both  filters  require  a  good  initial  estimate  of  the  float  positions  and  velocities  in  order 
perform  well.  Their  respective  responses  to  bad  initial  estimates,  which  are  discussed  below, 
include  becoming  stuck  at  a  bad  position  estimate,  slowly  converging  to  a  good  position  esti¬ 
mate,  and  no  effect  at  all.  Nearly  perfect  initial  estimates  of  the  floats’  initial  positions  and 
velocities  were  u.sed  in  the  simulations  presented  in  Appendix  C  in  order  to  isolate  the  effects  of 
changes  in  process  noise,  measurement  noise,  and  deployment  geometry.  In  this  Section,  a 
method  is  developed  for  obtaining  sufficiently  good  initial  estimates  of  float  positions  and  veloci¬ 
ties  to  enable  both  filters  to  produce  good  position  estimates  beginning  almost  immediately. 

The  filter’s  perform  differently  when  given  poor  initial  estimates,  A  really  bad  initial 
estimate  can  cause  the  least  squares  filter  to  converge  to,  and  remain  stuck  at,  the  wrong 
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Figure  4.1  Coordinate  system  for  the  Swallow  float  localization 


position.  This  position  is  called  a  local  mtntmum  because  it  is  a  local  minimum  in  the  loss  func¬ 
tion,  but  not  the  global  minimum  (see  Equation  3.9  of  Appendix  C).  However,  if  we  provide  the 
least  squares  filter  with  only  a  moderately  poor  initial  estimate,  so  that  it  does  not  get  stuck  at 
a  local  minimum,  then  its  first  position  estimate  will  be  very  good.  In  terms  of  the  estimate 
error,  a  really  bad  initial  estimate  can  result  in  a  large  initial  error  which  remains  large  for  sub¬ 
sequent  measurements.  However,  a  moderately  poor  initial  estimate  will  result  in  a  small  initial 
error  which  remains  small. 

The  Kalman  filter’s  response  to  a  poor  initial  estimate  depends  upon  two  other  param¬ 
eters,  the  estimated  initial  position  error  and  the  estimated  initial  velocity  error,  which  are  sup¬ 
plied  to  the  filter  separately.  If  the  estimated  initial  position  error  is  large,  then  the  filter  relies 
primarily  on  the  measurement  to  estimate  float  positions.  Likewise  for  the  initial  velocity  error. 
If  the  measurement  is  noisy,  the  filter’s  estimate  may  jump  around  considerably.  The  error  in 
the  filter’s  estimate  depends  entirely  upon  the  measurement  error. 

If  the  initial  position  and  velocity  error  estimates  are  set  small,  the  Kalman  filter  will 
estimate  float  positions  and  velocities  b.ised  upon  the  process  and  measurement  noise  estimates 
which  it  has  been  given.  If  there  was,  in  fact,  error  in  the  initial  estimate,  the  filter  will  con¬ 
verge  very  slowly  to  the  correct  position  estimate.  In  terms  of  the  estimate  error,  a  bad  initial 
estimate  causes  a  large  initial  error  which  slowly  decreases.  On  the  other  hand,  a  good  initial 
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estimate  results  in  a  small  initial  error  which  stays  small. 

To  summarize,  the  least  squares  filter  may  respond  to  a  very  bad  initial  estimate  by 
becoming  stuck  at  a  local  minimum.  However,  it  will  converge  immediately  to  a  good  position 
estimate  if  the  initial  estimate  is  not  too  bad.  The  Kalman  filter’s  response  to  a  poor  initial 
estimate  depends  upon  the  values  of  the  estimated  initial  position  and  velocity  errors.  If  they 
are  large,  the  filter  depends  entirely  upon  the  measurements.  If  they  are  small,  the  filter  will 
slowly  converge  to  the  correct  position  estimate. 

Good  initial  estimates  can  be  obtained  by  taking  advantage  of  the  strengths  of  both 
filters.  The  following  procedure  was  found  to  work  well  in  simulations; 

(1)  Starting  with  arbitrary  initial  float  position  and  velocity  estimates  and  large 
values  for  the  estimated  initial  position  and  velocity  errors,  the  Kalman  filter 
is  run  on  the  first  2  to  4  measurements.  Using  the  filter’s  estimates  after  the 
last  measurement  as  the  new  initial  estimates,  the  filter  is  again  run  on  the 
same  set  of  measurements.  This  step  is  repeated,  each  time  using  the  filter’s 
estimate  at  the  end  of  the  previous  run  as  the  starting  point  for  the  next  run, 
until  the  filter’s  position  estimate  stops  changing  much. 

(2)  Using  the  position  estimate  obtained  in  (1)  above,  the  least  squares  filter  is 
run  on  the  first  measurement.  This  produces  the  best  estimate  of  float  posi¬ 
tions  for  that  measurement  and  becomes  the  initial  position  estimate  for  both 
filters. 

(3)  Using  the  result  of  (2)  as  the  initial  position  estimate,  an  arbitrary  initial 
velocity  estimate,  a  small  estimated  initial  position  error  and  a  large 
estimated  initial  velocity  error,  the  Kalman  filter  is  again  run  on  the  first  2  to 
4  measurements.  The  filter’s  position  estimate  does  not  change  much  but  the 
velocity  estimate  does.  Using  the  filter’s  velocity  estimate  after  the  last  meas¬ 
urement,  the  filter  is  again  run  on  the  same  set  of  measurements.  This  step  is 
repeated,  each  time  using  the  filter’s  velocity  estimate  at  the  end  of  the  previ¬ 
ous  run  as  the  starting  point  for  the  next  run,  until  the  filter’s  velocity  esti¬ 
mate  stops  changing  much. 


4.1.3.  Relationship  Between  Actual  and  Predicted  Error  and  the  Innovations  Sequence 

The  performances  of  the  Kalman  and  least  squares  filters  are  compared  in  Appendix  C 
on  the  basis  of  the  rms  error  in  their  position  estimates.  However,  both  filters  also  estimate  the 
error  in  their  output.  We  might  ask  the  question  "What  is  the  relationship  between  the  error 
estimates  produced  by  the  filters  and  the  actual  position  errors?"  The  answer  to  this  question 
will  permit  us  to  predict  the  estimate  error  with  some  confidence  when  we  process  experimental 
data. 

The  least  squares  filter  produces  an  estimate  of  the  error  (actually  the  mean  magnitude 
of  the  residual),  which  is  defined  in  equation  (3.25)  of  Appendix  C.  It  is  analogous  to  the  mean 
square  difference  between  the  measurement  and  position  estimate,  normalized  by  the  number  of 
measurements. 

The  Kalman  filter  produces  an  estimate  of  the  position  error  covariance,  identified  in 
equation  (4.29)  of  Appendix  C.  Also  produced  is  the  mean  magnitude  of  the  innovation,  which  is 
a  measure  of  the  difference  between  the  predicted  float  positions  and  the  positions  indicated  by 
the  measurement,  and  is  identified  below  equation  (4.13)  in  Appendix  C.  The  latter  is  similar  to 
the  residual  produced  by  the  least  squares  filter. 

The  error  estimates  produced  by  the  filters  are  now  compared  to  the  actual  rms  errors 
for  the  simulations  presented  in  Appendix  C.  Figure  4.2  compares  the  predicted  and  actual  rms 
error  in  the  least  squares  filter  result  for  the  7  Cases  simulated.  The  filter’s  estimate  of  the 
error  comes  quite  close  to  the  actual  error  in  Cases  1,  2  and  7,  where  the  deployment  geometry 
was  normal  and  either  or  both  the  process  and  measurement  noise  were  low.  The  actual  and 
estimated  errors  are  higher  in  Cases  2  and  6  than  in  Cases  1  and  7  because  the  measurement 
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noise  is  higher  in  Cases  2  and  6.  The  filter  slightly  underestimates  the  error  in  Case  6,  when  the 
process  and  measurement  noise  are  higher.  Cases  3,  4,  and  5  show  that  the  filter  underestimates 
the  error  when  the  deployment  geometry  becomes  less  optimal.  Under  these  circumstances  the 
error  estimate  is  about  a  factor  of  2  lower  than  the  actual  error. 
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Figure  4.2  Predicted  and  actual  rms  error  in  the  least  squares  filter  results 
for  the  seven  simulated  cases 


Figures  4.3  and  4.4  compare  the  mean  magnitude  of  the  innovation  and  the  predicted 
and  actual  rms  error  obtained  using  the  Kalman  filter  in  the  4  Cases  involving  changes  in  pro¬ 
cess  and  measurement  noise.  The  vertical  axes  are  rms  error  or  mean  innovation  magnitude,  in 
meters.  The  horizontal  axes  are  the  estimate  of  process  noise  given  to  the  filter  as  a  parameter. 
The  process  noise  estimate  is  used  to  control  the  filter.  It  is  denoted  q  and  will  be  the  subject  of 
considerable  discussion  below. 

Each  panel  in  Figures  4.3  and  4.4  contains  3  curves.  The  top  curve  is  the  mean  inno¬ 
vation  magnitude  and  is  indicated  by  |  /  |  .  The  middle  curve  is  the  true  rms  error,  and  is 
labeled  e .  The  bottom  curve  is  the  estimate  of  the  standard  deviation  of  the  error  produced  by 
the  filter  and  is  labeled  e. 

The  Cases  shown  in  Figure  4.3  have  low  underlying  process  noise.  The  actual  error  is 
bracketed  by  the  innovation  and  the  predicted  error,  and  their  relative  values  are  approxi¬ 
mately  constant.  The  error  in  the  estimate  varies  with  q:  as  q  is  lowered,  and  the  error 
decreases,  and  vice  versa.  The  lowest  error  is  obtained  by  setting  q  very  low. 

Increasing  the  measurement  noise  increases  the  spacing  between  the  3  curves,  as  may 
be  seen  by  comparing  the  upper  and  lower  plots  in  Figure  4.3. 

The  Cases  shown  in  Figure  4.4  show  the  innovation  and  rms  error  when  the  underlying 
process  noise  is  high.  When  q  is  set  large,  the  curves  resemble  the  low  process  noise  cases. 
When  q  is  set  low,  however,  the  innovation  and  actual  error  rise  but  the  estimated  error 
remains  low.  This  phenomenon  is  called  divergence  because  the  filter’s  estimate  of  position  error 
diverges  from  the  true  position  error. 
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Figure  4.4  Mean  magnitude  of  the  innovation,  predicted  and  actual  rms 
error  using  the  Kalman  filter  on  simulation  Cases  6  and  7 
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Figures  4.3  and  4.4  show  that  the  error  in  the  Kalman  filter  estimate,  the  mean  inno¬ 
vation  magnitude  and  the  estimated  error  vary  as  a  function  of  q,  the  estimate  of  process  noise 
given  to  the  filter.  Thus  q  is  the  parameter  which  is  used  to  control  the  filter’s  performance. 
The  filter  uses  q  to  determine  how  much  weight  to  give  to  its  own  track  of  the  float  positions, 
relative  to  the  measurements.  When  q  is  set  high,  indicating  that  the  process  noise  is  thought 
to  be  high,  the  filter  gives  greater  weight  to  the  measurements  and  less  to  its  own  track  of  the 
floats.  In  this  case,  the  Kalman  filter  result  approaches  that  of  the  least  squares  filter  in  that 
the  position  error  depends  upon  the  measurement  noise. 

When  q  is  set  low,  indicating  that  the  process  noise  is  thought  to  be  low,  the  filter 
gives  greater  weight  to  its  own  track  of  the  floats  and  less  to  the  new  measurements  coming  in. 
If  the  process  noise  really  is  low,  the  filter  works  fine.  If  the  process  noise  is  actually  high,  how¬ 
ever,  divergence  can  occur  and  the  po.sition  error  becomes  large.  This  is  seen  when  q  <  lO'^  in 
Figure  4.4. 

The  optimal  value  of  q,  which  minimizes  the  rms  error  in  the  filter’s  estimate,  is  readily 
apparent  in  the  simulations  shown  in  Figures  4.3  and  4.4  because  the  true  rms  position  error  is 
known.  Figure  4.4,  however,  shows  how  to  determine  q  when  processing  experiment  data  and 
the  true  position  error  in  unknown.  The  optimum  value  for  q  is  seen  to  correspond  to  the 
minimum  mean  innovation  magnitude.  This  suggests  that  when  experimental  data  are  being 
processed,  the  filter  can  be  run  several  times  using  different  values  for  q,  and  the  optimum  value 
of  q  selected  as  the  one  which  minimizes  the  innovation. 

The  error  in  the  resulting  position  estimate  will  be  bracketed  by  the  mean  magnitude 
of  the  innovation  and  the  estimated  standard  deviation  of  the  error.  A  reasonable  estimate  for 
the  rms  position  error  would  be  the  average  of  the  mean  innovation  and  estimate  of  rms  error. 
A  conservative  estimate  for  the  rms  error  would  be  the  mean  magnitude  of  the  innovation. 


4.1.4.  Effect  of  Sound  Speed  Errors  on  Location  Error 

The  simulations  in  Appendix  C  included  travel  time  error  but  no  sound  speed  error. 
What  remains  is  to  (1)  estimate  the  sound  speed  error  and  (2)  determine  its  impact  on  the 
filters’  performance. 


4. 1.4.1.  Sound  Speed  Errors 

The  sound  speed  estimate  is  usually  based  upon  ocean  temperature  versus  depth  meas¬ 
urements  made  at  the  experiment  site  using  expendable  bathythermographs  (XBT’s).  The  sound 
speed  profile  is  then  calculated  using  historical  values  for  salinity.  There  will  be  errors  in  the 
sound  speed  profile  and  it  will  fluctuate  somewhat  during  the  experiment.  Typical  measure¬ 
ments  of  maximum  sound  speed  fluctuations  reported  in  the  literature  are  0.75  and  1.0 
m/sec. 


An  estimate  of  the  errors  and  variation  found  in  typical  sound  speed  estimates  was 
obtained  using  a  series  of  12  XBT  measurements  made  during  a  26  hour  period  by  MPL  person¬ 
nel  from  FLIP  during  September  1987  at  35  *  N,  126  ’  W,  about  265  miles  west  of  Point  Concep¬ 
tion,  California.  The  standard  deviation  of  the  sound  speed  at  a  particular  depth  was  just 
under  0.6  meters/ sec.  This  number  is  taken  as  representative  of  the  rms  sound  speed  error  at  a 
particular  depth  for  experiments  lasting  approximately  24  hours. 

It  is  shown  in  Appendix  C  that  the  filters  utilize  the  average  sound  speeds  over  the 
paths  separating  each  pair  of  floats.  To  calculate  the  average  sound  speeds,  the  sound  speed 
measurements  are  first  used  to  develop  a  multi-layer  ocean  model  with  constant  sound  speed  in 
each  layer.  Then,  for  each  float  pair,  the  average  sound  speed  is  found  from 
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where  rf,  and  c,  are  the  thickness  and  sound  speed  for  the  layer  and  the  summation  is  over 
all  of  the  layers  between  the  float  pair. 

Equation  (4.2)  is  a  good  approximation  so  long  as  the  the  ray  paths  are  approximately 
straight.  The  Generic  Sonar  Model  (GSM)  ray  tracing  program^^  was  employed  to  determine  if 
that  condition  is  satisfied  during  the  September  1986  deployment.  Interfloat  travel  times  were 
calculated  for  representative  pairs  of  floats  using  the  sound  speed  profile  shown  in  Figure  3.12 
and  the  approximate  positions  of  the  floats  during  the  September  1986  experiment.  The 
effective  sound  speed  was  computed  by  dividing  the  geometric  distances  separating  the  floats  by 
the  GSM-produced  travel  times.  The  average  and  effective  sound  speeds  differed  by  less  that  0.5 
meters/sec,  indicating  that  equation  (4.2)  is  a  good  approximation. 

The  rms  error  of  the  average  sound  speed  c  is  a  complicated  function  of  the  errors  in 
the  sound  speeds  c,’s.  However,  since  (4.2)  is  essentially  an  average,  and  since  the  sound  speed 
errors  are  taken  to  be  approximately  mean  zero,  the  error  in  the  average  sound  speed  would  be 
expected  to  be  less  than  that  in  the  sound  speed  for  a  particular  layer. 


4.I.4.2.  Efifect  On  Position  Error 

In  Appendix  C,  the  measurement  noise  is  defined  to  be  the  error  in  the  interfloat  or 
float  to  surface  ranges,  that  is,  to  include  the  error  made  in  measuring  the  travel  time  as  well  as 
the  error  made  in  estimating  the  sound  speed  along  the  path  over  which  the  range  pulse  travels. 
From  equation  (3.8)  of  Appendix  C,  the  rms  sound  speed  and  measurement  errors  for  floats  i 
and  i  are  related  according  to 


(4.3) 


where  =  0.6  m/sec  is  the  estimated  rms  error  in  the  average  sound  speed,  c  is  the  average 
sound  speed  estimate,  f  is  the  travel  time  estimate,  a,  is  the  rms  travel  time  error  estimate,  and 
(T,  is  the  rms  measurement  error.  It  is  shown  in  (21)  that  the  rms  travel  time  error  during  the 
September  1986  experiment  is  approximately  0.0013  seconds.  Take  the  average  sound  speed  to 
be  1500  m/sec.  If  the  path  length  between  floats  i  and  j  is  approximately  3  km,  so  that  the 
travel  time  is  about  2  seconds,  then  the  rms  measurement  error  is 


r  l'/2 

=  1(4  sec^)  (0.36  m^sec^)  +  (2.25  X  10*  mVsec*)  (1.69  X  lO"^  sec^)  I 

1.44  -»-3.8j*/^  meters  =2.3  meters  (r  =  2  seconds)  (4.4) 

The  second  term,  which  involves  the  rms  travel  time  error,  is  about  twice  as  large.  If  the  range 
is  increased  to  7.5  km,  so  that  the  travel  time  increases  to  5  seconds,  then 

<7^  =  ^9  +  3.8  meters  =  3.6  meters  (r  =  5  seconds) 

The  first  term  now  dominates  and  the  rnis  measurement  error  is  increased  by  about  2/3. 

The  effect  of  sound  speed  errors  is  to  increase  measurement  noise,  which  we  have  seen 
increases  error  in  the  filters’  results.  However,  the  increase  in  measurement  noise  is  apparently 
not  significant  except  for  the  largest  travel  times.  Since  both  filters  accept  estimates  of  the  rms 
sound  speed  and  travel  time  error  and  use  equation  (3.8)  of  Appendix  C  to  calculate  measure¬ 
ment  error  variance,  they  gives  less  weight  to  longer  ranges  and  thus  properly  respond  to  sound 
speed  errors. 


(4.5) 
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4.2.  Application  to  Experiment  Data 

September  1986  experiment  acoustic  range  data  between  records  1150  and  1870  are 
now  applied  to  the  least  squares  and  Kalman  filters.  A  coordinate  system  is  established  in 
which  the  origin  lies  at  the  surface  directly  above  float  0,  and  float  1  lies  in  the  XZ  plane.  This 
is  shown  in  Figure  4.5.  The  Z  axis  is  vertical  and  extends  through  float  0.  The  directions  of  the 
X  and  Y  axes  relative  to  magnetic  north  are  determined  from  Figure  3.3  to  be  approximately 
280“  and  190°,  respectively. 


Magnetic 

North 


Figure  4.5  Plan  view  of  the  coordinate  system  established  for  localizing 
Swallow  floats  during  the  September  1986  experiment 


The  processing  steps  shown  in  Figure  4  of  Appendix  C  were  followed  to  obtain  travel 
times  estimates  between  each  float  pair  and  between  each  float  and  the  surface.  The  rms  travel 
time  error  was  also  estimated  for  each  pair  and  for  the  surface  echoes.  The  average  rms  error  is 
found  to  be  1.3  msec. 

The  procedure  described  in  Section  4.1.1  was  used  to  obtain  initial  position  and  velo¬ 
city  estimates.  Using  an  rms  sound  speed  error  of  0.6  m/sec,  the  travel  times  were  applied  to 
the  least  squares  and  Kalman  filters,  and  the  results  shown  in  the  upper  panel  of  Figure  4.6. 
The  mean  innovation  magnitude  and  predicted  standard  error  produced  by  the  Kalman  filter 
are  labeled  |  I  |  and  e*,  respectively.  The  mean  error  predicted  by  the  least  squares  filter  is 
indicated  by  a  dotted  line  and  labeled  Cj,. 

The  innovation  magnitude  and  predicted  rms  error  look  much  like  the  high  process 
noise  plots  in  Figure  4.4,  except  shifted  up  and  to  the  right.  The  right  shift  suggests  that  the 
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Figure  4.6  Results  of  the  Kalman  and  least  squares  filters  applied  to  the 
September  1986  experiment  an  simulation  Case  1 1 
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process  noise  in  the  experiment  data  was  higher  than  that  used  in  the  simulations,  and  the  shift 
up  indicates  likewise  for  th“  measurement  noise.  However  the  rms  travel  time  error  found  in  the 
experiment  data  (1.3  msec)  is  approximately  equal  to  that  used  in  the  high  measurement  noise 
simulations  (1.4  msec).  The  shift  up,  then,  is  apparently  due  to  sound  speed  errors  in  the  experi¬ 
mental  data. 

In  order  to  match  the  experimental  results  more  closely,  another  simulation  was  run  in 
which  the  rms  process  noise  was  increased  by  a  factor  of  about  30  and  the  rms  travel  time  error 
increased  from  about  1.4  to  1.7  msec.  The  results  are  shown  in  the  bottom  panel  of  Figure  4.6, 
labeled  as  in  the  upper  panel  except  that  e*  and  Cj,  are  the  true  rms  errors  in  the  estimates. 
The  Case  11  Kalman  filter  results  are  quite  similar  to  the  experimental  results.  The  mean  inno¬ 
vation  and  predicted  error  have  approximately  the  same  shape  and  location  in  both  plots.  The 
optimal  value  of  q  appears  to  be  10“*  for  the  experiment  data  and  10“*°  for  Case  11.  The 
difference  between  these  two  numbers  could  be  a  result  of  the  order  of  the  simulation  model 
being  too  low. 

Based  upon  the  mean  innovation  magnitude  and  predicted  standard  deviation  of  the 
error,  the  ims  error  in  the  Kalman  filter’s  position  estimate  for  the  September  1986  experiment 
is  estimated  to  be  about  2.4  meters.  It  is  certainly  less  than  3.5  meters.  The  rms  difference 
between  the  Kalman  filter  and  the  least  squares  position  estimates  was  found  to  be  about  3 
meters.  Based  upon  the  error  estimates  produced  by  both  filters,  and  given  the  apparent  high 
process  noise  in  the  experimental  data,  it  is  impossible  to  determine  whether  one  estimate  is 
better  than  the  other.  Both  apparently  satisfy  the  condition  in  (4.1).  The  Kalman  filter  posi¬ 
tions  were  chosen  for  the  beamforming  in  Chapter  5. 

Plan  views  of  the  float  positions  as  a  function  of  time  are  shown  in  Figures  4.7  and  4.8. 
The  former  shows  all  of  the  floats  and  the  latter  shows  only  the  freely-drifting  floats.  It  is 
interesting  to  note  that  four  of  the  floats  (3,  4,  5  and  10)  drifted  in  a  generally  counter-clockwise 
direction,  while  float  9  drifted  generally  in  a  generally  clockwise  direction.  According  to  Table 
3.1,  the  four  floats  which  drifted  together  all  deployed  to  about  1000  meters  depth,  while  float  9 
deployed  to  about  800  meters  depth. 

All  five  of  the  floats  drifted  with  an  average  speed  of  approximately  200  meters/ hour, 
which  is  about  0.1  kt  or  5  cm/sec. 


5.  Beamforming  With  the  Swallow  Float  Array 

Having  determined  the  float  positions  during  the  September  1986  experiment,  the  mid¬ 
water  column  floats  may  now  be  treated  as  an  array  of  sensors  and  their  time  series  beam- 
formed.  In  general,  arrays  of  spatially  separated  sensors  are  employed  in  measuring  a  sound 
field  in  order  to  determine  its  directionality,  which  is  the  direction  from  which  sound  is  arriving. 
This  is  the  goal  of  beamforming  with  the  Swallow  floats. 

Beamforming  with  an  array  of  Swallow  floats  turns  out  to  be  considerably  different 
from  beamforming  with  an  array  of  omni-directional  hydrophones.  This  is  because  each  float 
contains  3  orthogonally-oriented  geophones,  and  each  geophone’s  response  depends  upon  the 
angle  between  the  source  direction  and  the  geophone  axis.  Simulations  are  used  to  demonstrate 
differences  between  the  results  obtained  using  omni-directional  and  directional  sensors.  The 
experimental  results  are  thus  interpretable,  although  the  large  array  aperture  and  an  inhomo¬ 
geneous  sound  field  cause  them  to  be  somewhat  confusing. 
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Figure  4.7  Kalman  filter  estimate  of  time-varying  float  positions  between  records  11  SO  and  1870 
(approximately  9  hours  duration)  fcH’  the  September  1986  experiment.  Floats  0, 1  and  2  are 
on  the  bottom;  floats  3, 4, 5. 9  and  10  are  freely-drifting.  "O"  marks  the  starting  position. 
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Figure  4.8  Kalman  filter  estimate  of  time-varying  float  positions  between  records  1150  and  1870 
(apiHDximately  9  hours  duration)  for  the  Sq)tember  1986  experiment.  Dq}ths  are  aiq)ioxi- 
mately  1000  meters  for  floats  3, 4, 5,  and  10,  and  approximately  800  meters  for  float  9.  "O" 
marlcs  the  starting  position.  ”+”  marks  3  hour  intervals. 
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5.1.  Theoretical  Development  and  Simulation  Results 

The  job  of  a  time  delay  and  sum  beamformer  is  to  combine  array  element  signals  to 
produce  beams  which  listen  preferentially  in  particular  directions,  called  look  directions.  For 
each  look  direction,  a  set  of  delays  must  be  calculated  which,  when  applied  to  the  element  time 
series,  cause  signals  arriving  from  that  direction  to  be  aligned  in  phase,  or  co-phased.  The 
amplitude  of  the  sum  of  harmonic  signals  with  the  same  phase  is  the  sum  of  their  amplitudes. 
This  is  greater  than  the  sum  of  harmonic  signals  which  do  not  have  the  same  phase. 

Consider  an  array  of  sensors  and  one  or  more  harmonic  sources.  Assume  that  all  of  the 
sources  are  located  far  from  the  array  so  that  wavefront  curvature  is  minimal  over  the  array 
aperture.  As  defined  in  Figure  5.1,  the  azimuthal  bearing  from  the  array  to  the  source  is  meas¬ 
ured  counterclockwise  from  the  X  axis  and  denoted  by  9,.  The  source  elevation  angle  is  meas¬ 
ured  from  the  horizontal  and  denoted  If  the  array  is  excited  by  a  complex  exponential  sig¬ 
nal 

«(0  =  e‘“‘  (51) 

whose  bearing  is  9,,  4>,  from  the  array,  then  the  beamformer  output  can  be  written 

b{t)  =bT(9i,4>u  9 „  <!>.)■  (5.2) 

where  the  array  bearing  response  bri^i,  <j>i,  9,,  d>,)  is  a  complex  quantity  which  contains  the 
amplitude  and  phase  of  the  response."®  This  is  a  useful  notation  because  the  bearing  responses 
of  different  sensors  and  arrays  can  be  compared  to  each  other  without  reference  to  the 
term. 


5.1.1.  Single  Element  Arrays 

Consider  first  a  simple  array  composed  of  a  single,  omni-directional  element.  For  sim¬ 
plicity,  start  with  a  source  which  is  in  the  XY  plane  (<l>,  =0)  and  constrain  the  look  direction  to 
also  lie  in  the  XY  plane  (<t>i  =0).  Since  the  array  contains  only  one  sensor  whose  response  does 
not  depend  upon  the  direction  to  the  source,  the  array  bearing  response  does  not  depend  upon 
either  the  source  or  look  direction  azimuth  angles. 

br{9i,  6,)  =  1  (single  omni-directional  element)  (5.3) 

Now  substitute  a  single  geophone  for  the  omni-directional  sensor.  The  geophone’s 
response  is  directional,  i.e.  it  depends  upon  the  angle  between  its  axis  and  the  bearing  to  the 
source.  Denote  the  geophone  azimuth  angle  relative  to  the  array  X  axis  as  9/f  as  shown  in  Fig¬ 
ure  5.2,  and  set  the  geophone  elevation  angle  4>h  equal  to  zero.  The  geophone  bearing  response 
is 

br{0i,  9,)  =cos  {9if  —9,)  (single  geophone)  (5.4) 

This  is  called  a  dipole  response  because  it  has  the  same  form  as  the  field  of  an  electric  dipole. 

If  the  geophone  axis  is  rotated  until  it  is  parallel  to  the  X  axis  of  the  array,  its  bearing  response 
becomes 

br(9i,  0,)  =cos  (5,)  (single  geophone  parallel  to  X  axis)  (5.5) 
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Figure  5.1  Coordinate  system  used  for  beamforming 
0  -  azimuth  angle  in  the  XY  plane 
0  -  elevation  angle  out  of  the  XY  plane 


Each  Swallow  float  contains  three  orthogonally-oriented,  collocated  geophones  whose 
axes  are  assumed  for  now  to  coincide  with  the  X,  Y  and  Z  axes  of  the  array.  For  the  geophone 
whose  axis  is  parallel  to  the  Y  axis,  0f{  =90  '  and  its  bearing  response  is; 


®.)  =cos  (90  —0,)  =sin  (0,)  (single  geophone  parallel  to  Y  axis)  (5.6) 


Now  consider  an  array  of  two,  collocated  geophones,  one  parallel  to  the  X  and  the 
other  parallel  to  the  Y  axis,  as  shown  in  Figure  5.3.  The  array  bearing  response  can  be  maxim¬ 
ized  in  the  direction  of  the  source  by  forming  a  time  series  for  each  look  direction  in  which  each 
sensor  signal  is  multiplied  by  the  cosine  of  the  angle  between  its  axis  and  the  look  direction  (the 
look  direction  cosines). 
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Substituting  (5.5)  and  (5.6)  for  the  bearing  responses  of  the  individual  geophones  and  using  cos 
(90 '  -  Of)  =  sin  (9() 


=  cos  {9,  —9,)  (X  and  Y  axis  geophones)  (5.7) 


It  is  clear  from  equation  (5.7)  that  the  bearing  response  of  the  collocated,  orthogonally-oriented 
geophone  pair  can  provide  information  about  the  source  azimuth,  whereas  the  omni-directional 
element  bearing  response  of  (5.3)  does  not.  Comparing  the  bearing  response  of  two  geophones  to 
that  of  one  omni-directional  sensor  may  seem  unfair.  However,  the  geophone  pair  are  packaged 
as  one  unit,  so  that  it  is  actually  the  bearing  responses  of  two  sensor  packages  which  are  being 
compared. 

This  is  an  important  comparison  because  a  major  part  of  the  effort  in  measuring  sound 
in  the  ocean  is  in  getting  the  sensor  in  place,  getting  power  to  it,  storing  its  output,  and  retriev¬ 
ing  it.  Deploying  two  separate  packages  requires  very  nearly  twice  the  hardware  and  effort  of 
deploying  one  package.  Deploying  two  geophones  packaged  as  one  unit  requires  about  the  same 


37 


R.  L.  Culver 


Y 


Figure  6.3  Geometry  for  two  collocated  geophones,  one  parallel  to  the 
X  axis  and  one  parallel  to  the  Y  axis 


hardware  and  effort  as  deploying  one  omni-directional  element.  Thus  the  advantage  accrued  by 
deploying  geophone  pairs  (or  as  we  will  see  shortly,  geophone  triplets),  is  that  source  bearing 
information  is  obtained  with  just  a  single  package. 

If  there  is  more  than  one  source  transmitting  a  signal  to  the  array,  then  the  bearing 
response  of  the  single  omni-directional  sensor  is  the  sum  of  the  responses  to  the  sources  individu¬ 
ally,  which  is  independent  of  the  source  bearing  angles. 

®»i)  +  ®«2)  =  2  (single  omni-directional  element)  (5.8) 

Similarly,  the  bearing  response  of  the  collocated,  orthogonally-oriented  geophone  pair  is  the  sum 
of  the  responses  to  the  sources  individually: 

br(&t,  ^,,)  -I-  br($i,  0,2)  =  cos  (Sj  -  6,^)  +  cos  (5j  —  0,2)  (X  and  Y  axis  geophones)  (5.9) 

Equation  (5.9)  has  a  peak  when  the  look  direction  equals  either  of  the  source  directions. 
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5.1.2.  Multi-Element  Arrays 

Now  consider  an  array  of  sensors.  The  elements  are  uniformly  spaced  on  the  Y  axis  at 
100  meters  intervals  as  shown  in  Figure  5.4.  The  array  bearing  response  is  no  longer  indepen¬ 
dent  of  the  source  and  look  directions  because  the  elements  are  spatially  separated  and  the 
beamformer  can  co-phase  signals  coming  from  the  look  direction.  For  simplicity  in  the  develop¬ 
ment,  the  source  and  look  direction  continue  to  be  constrained  to  the  5^  plane.  The  results 
will  be  expanded  to  three  dimensions  subsequently. 


Y 


Figure  5.4  Multi-element  array  geometry 


Figure  5.5  compares  the  magnitude  of  the  bearing  responses  of  the  omni-directional  ele¬ 
ment  and  orthogonally-oriented  geophone  arrays  with  a  6.7  Hz  source  at  =0  ".  Both  bearing 
responses  have  peaks  in  the  source  direction,  which  is  desirable.  However,  both  bearing 
responses  also  have  peaks  in  another  direction  too,  which  indicates  directional  ambiguity  with 
respect  to  the  array  axis  (the  Y  axis  of  the  coordinate  system).  When  the  source  bearing  is  0  ° 
(source  broadside  to  the  array),  we  know  that  it  is  broadside  but  we  cannot  determine  which 
side  it  is  on  (0°  or  180° ).  This  is  known  as  the  right-left  ambiguity  of  a  line  array  because  if  we 
stand  at  the  origin  and  look  down  the  array  axis,  we  cannot  tell  rather  the  source  is  coming 
from  our  right  or  left.  The  ambiguity  (0°  or  180°)  persists  in  the  geophone  array  bearing 
response  because  the  magnitude  of  the  geophone’s  cos  (fl;)  response  has  a  180”  (forward- 
backward)  ambiguity,  which,  when  the  source  angle  is  0°,  coincides  with  the  line  array’s  right- 
left  ambiguity. 
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Beorlng  Resporpe  of  4  OiiinL-Djrectlcnol  Elements  on  the  Y  Axis 
b^lr^  lE  100  nmeTEi  o.7  Hz  eorce  ol  ffelmuth  0  Degrees 


Bearing  Respprse  of  4  Geoghqne  Pofrs  on  the  Y  fixl9„  ^ 
Spacing  is  100  raetersi  6,7  Hz  source  ot  Azimuth  60  Degrees 


Figure  5.5  Comparison  of  the  magnitude  of  the  bearing  rc-sponscs  of  omni-directional  clement  and  geo- 
phone  arrays  with  a  6.7  Hz  source  at  0,  =  (Y’.  The  array  elements  ate  uniformly  spaced  at 
100  meter  intervals. 
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The  half-wavelength  spacing  for  this  frequency,  with  c  =  1500  meters/sec,  is 

X  ^  ^  ^  1500  meters/sec  ^  (5.10) 

2  2  2  6.7  Uz  '  ' 

Thus  the  element  spacing  of  100  meters  used  in  the  simulations  is  less  than  X/2.  This  is  why  in 
the  omni-directional  bearing  response,  the  two  side  lobes  merge  at  0;  =  90  “  and  6i  =  270  ” 
rather  than  going  to  zero  as  they  would  if  the  element  spacing  were  exactly  X/2. 

Note  from  Figure  5.5  that  the  geophone  array  achieves  about  5  dB  side  lobe  level 
reduction.  This  is  because  of  the  directionality  of  the  geophones  parallel  to  the  X  axis,  and 
because  the  geophones  parallel  to  the  Y  axis  do  not  detect  the  signal  at  all. 

Figure  5.6  compares  the  magnitude  of  the  bearing  responses  of  the  omni-directional  ele¬ 
ment  and  orthogonally-oriented  geophone  arrays  with  a  6.7  Hz  source  at  S,  =60  '.  The  direc¬ 
tional  ambiguity  in  the  omni-directional  element  bearing  response  is  such  that  a  source  bearing 
60“  from  the  array  cannot  be  distinguished  from  one  bearing  120“ .  This  ambiguity  is  reduced 
by  5  dB  in  the  geophone  array  bearing  response  because  the  geophone  and  array  ambiguities  are 
not  pointed  in  the  same  direction.  The  geophone  array  is  again  seen  to  achieve  considerable 
side  lobe  reduction. 


5.1.3.  Sparse  Multi-Element  Arrays 

The  advantages  offered  by  Swallow  floats  in  closely-spaced,  uniform  arrays  are  clear 
enough,  but  we  saw  in  Chapter  4  that  the  mid-water  column  floats  were  1000  to  2200  meters 
apart  during  the  September  1986  experiment.  Arrays  whose  element  spacing  is  greater  than 
one-half  the  wavelength  of  sound  at  the  highest  frequency  of  interest  are  said  to  be  sparse.  The 
bearing  response  of  sparse,  uniformly-spaced  arrays  may  contain  lobes  whose  levels  are  as  high 
as  main  lobe,  in  directions  other  than  the  source  direction.  They  are  called  grating  lohes!^ 

Figures  5.7  and  5.8  compare  the  bearing  responses  of  omni-directional  element  and  geo¬ 
phone  arrays  with  a  350  meter  element  spacing  and  with  6.7  Hz  sources  at  0“  and  60".  The 
numerous  grating  lobes  make  it  impossible  to  determine  the  source  direction  using  the  omni¬ 
directional  element  array.  The  geophone  array  reduces  grating  lobes  and  makes  it  much  easier 
to  recognize  the  source  bearing.  The  180“  ambiguity  persists  with  the  source  at  0“  because  the 
X  axis  geophones  and  the  array  both  have  ambiguities  in  this  direction.  The  source  direction 
ambiguity  is  reduced  5  dB  when  the  source  is  at  60  “ . 


5.1.4.  Three-Dimensional  Beamforming 

The  comparison  made  thus  far  in  this  Chapter  is  now  expanded  to  include  look  direc¬ 
tions  and  source  bearings  not  restricted  to  the  XY  plane.  As  in  equation  (5.3),  the  single,  omni¬ 
directional  element  bearing  response  is  independent  of  look  or  source  bearing  angles 

br(0,,  <fii^  6,  0,)  =  1  (single  omni-directional  element)  (5. 11) 


Also  as  in  the  two-dimensional  case,  the  bearing  response  of  the  geophone  depends 
upon  the  angle  between  the  gcophone  axis  and  the  source  direction.  If  that  angle  is  denoted  a, 
then 


br{0i,  ibi  0,  0,)  =  cos  (o)  (single  omni-directional  element)  (5.12) 


The  definition  of  the  dot  product  can  be  used  to  obtain  an  expression  for  cos  (o)  in  terms  of  Oi, 
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Figure  5.6  Comparison  of  ihc  magnitude  of  the  bearing  responses  of  omni-dirccUonal  element  and  gco- 
phone  arrays  with  a  6.7  Hz  sourep  at  9,  =  60®.  The  array  elcmenu  arc  uniformly  qtaced  at 
100  meter  intervals. 
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Beoring  Resooree  of  4  OmnL-pirectlonol  on  ^  Y  Axis 

S^Ir^  IB  350  motersi  6.7  Hz  sarco  at  R  Imuth  0  Degrees 


Bearing  Resogn^  of  4  Geoahcne  Pairs  on.tte  Y  Axis  _ 
Spacing  le  350  metersi  6.7 Ttz  source  at  Azimuth  0  Degrees 


Figure  5.7  Comparison  of  the  magnitude  of  the  bearing  responses  of  omni-directional  clement  and  geo- 
phone  arrays  with  a  6.7  Hz  source  at  0,  =  0°.  The  array  elements  are  uniformly  spaced  at 
350  meter  intervals. 
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Figure  5.8  Comparison  of  the  magnitude  of  the  bearing  Fc.eponses  of  omni  -directional  element  and  gco- 
phonc  arrays  with  a  6.7  Hz  source  at  0,  s  60*'.  The  array  elements  arc  uniformly  qxaced  at 
350  meter  intervals. 
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0,,  and  \{  U  and  ?  are  unit  vectors  in  the  geophone  axis  and  source  directions,  respec¬ 
tively,  then 

^  =  I  ^  I  I  ?  I  cos  (a)  (5.13) 

Denoting  the  geophone  azimuth  and  elevation  angles  6n  and  <l>n,  respectively,  then  using  the 
coordinate  system  defined  in  Figure  5.1,  the  components  of  are 

/?  =cos  {0„)  cos(<^h)  ^  +sm  {0„)  cos(<l>„)  ej  +sin[(t>ff)  ^  (514) 

where  e^,  and  e7  are  unit  vectors  along  the  axes.  Similarly,  the  components  of  ?  are 

fj  =cos  (^5)  cos(<t>s)  -I- sin  (O5)  cos(^s)  ^  +sin(^5)  eX  (5.15) 

The  dot  product  is  thus 

cos  (a)  =  FI  ■  ^  —  cos  (Ojf)  cos{4>h)  cos  (O^)  cos(<i>s)  -t-sin  (Ofj)  cos[4>fj)  sin  (Og)  cos((t>s) 

+  sin(<j>n)  sin(<^s) 

so  that 

^«)  =COS  (Off  —d,)  cos{<l>n)  cos(^,)  +Sll\(<j>n)  sin(^,)  (single  geophone)  (5.16) 

If  the  geophone  is  rotated  until  its  axis  is  parallel  to  the  array  X  axis 
(9^  =  0  ‘ ,  <l)f/  =  0  '),  its  bearing  response  becomes 

br{9i,  01  <^,)  =cos  (0,)  cos  (0,) (single  geophone  parallel  to  X  axis)  (5.17a) 

Similarly,  the  bearing  response  of  a  geophone  parallel  to  the  Y  axis  {Off  =90°,  (t>ff  =0  °)  will  be 

br{9i,  0,  0,  (j),)  =sin  {0,)  cos  (<^,)  (single  geophone  parallel  to  Y  axis)  (5.17b) 

and  that  of  a  geophone  parallel  to  the  Z  axis  (^//  =  90  ”)  will  be 

br(0i,  4>i  9,  Ip,)  =sin  {<p,)  (single  geophone  parallel  to  Z  axis)  (5.17c) 

Now  consider  three  colioc^ed  geophones,  one  oriented  parallel  to  each  axis.  The  bear¬ 
ing  response  of  the  combination  is  formed  by  multiplying  the  sensor  signals  by  the  look  direction 
cosines,  which  are 

/,  =cos  (Of)  cos  (P,) 
ly  =sin  (9i)  cos  (ipi) 
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=siii  ((^.,)  (5.18) 

The  combined  bearing  response  is  thus 

'A/.  <^.)L  .  .  ■  cos  (0,)  cos  (0,) 

JX,  Y  and  Z  axis  geopboaes  L  JX  axis  geopnoDe 

+  ^  ■  sin  [9,)  cos  (4) i) 

JY  axis  geophone 

+  \br(9i,  4>iJ,,4>,)\  810(0,) 

L  JZ  axis  geophone 

=  cos  (9,)  cos  (0,)  cos  (0,)  cos  (0,)  +sin  (0,)  cos  (0,)  sin  ((?,)  cos  (0,)  +sin  (0,)  sin  (0,) 


=  cos  (9,  -  9i)  cos  (0,)  cos  (0,)  +  sin  (0,)  sin  (0,)  (5.19) 

Equation  (5.19)  has  a  maximum  when  the  source  and  look  directions  are  equal.  When  the 

source  and  look  direction  azimuth  angles  are  equal,  equation  (5.19)  reduces  to  cos  (0,  —  0,), 
which  has  a  maximum  when  the  look  and  source  direction  elevation  angles  are  equal. 

Figure  5,9  shows  the  bearing  response  of  a  four  omni-directional  element  array  to  a  6.7 
Hz  source  at  azimuth  9,  =  0  ‘  and  elevation  0,  =0  “.  The  array  elements  are  uniformly  spaced 
at  100  meter  intervals.  The  look  direction  azimuth  angle  ranges  from  0'  to  360*,  and  the 
elevation  angle  varies  from  -30*  to  30* .  Each  trace  indicates  the  bearing  response  for  a  partic¬ 
ular  look  direction  elevation  angle  0,.  The  vertical  scale  for  each  trace  is  -10  to  0  dB,  i.e.  the 
trace  is  clipped  at  -10  dB. 

Several  points  can  be  made  about  this  bearing  response.  First,  the  0,  =0  trace  is 
identical  to  the  top  10  dB  of  the  upper  plot  in  Figure  5.5,  which  was  also  made  with  0,  =0. 

Second,  the  levels  of  the  peaks  at  =  0  '  and  180  *  is  independent  of  0,.  This  is  due  to  the 

array’s  right-left  ambiguity,  which  becomes  conical  in  3  dimensions.  When  the  source  is  broad¬ 
side,  the  directional  ambiguity  is  circular  so  that  the  cone  angle  is  180° .  Third,  the  lobe  width 
increases  as  |  0,  |  increases.  (This  can  be  seen  by  looking  closely  at  the  Figure).  This  is  due 
to  the  coordinate  system  used  in  the  beamforming.  As  the  elevation  angle  increases,  and  as  we 
sweep  through  a  360  *  azimuth  angle,  an  increasingly  greater  percentage  of  the  sweep  is  within 
the  main  lobe.  When  the  look  elevation  angle  is  straight  up  or  straight  down  (0  =±  90  ",  not 
.shown  in  Figure  5.9),  the  entire  sweep  is  pointed  toward  the  peak  of  the  main  lobe  and  the  level 
of  the  trace  is  0  dB  independent  of  azimuth  angle. 

Now  replace  each  omni-directional  element  with  three  collocated,  orthogonally-oriented 
geophones,  one  parallel  to  each  axis.  Figure  5.10  shows  the  bearing  response  of  the  array  to  the 
same  source  as  in  Figure  5.9.  The  geophone  array  response  appears  quite  similar  to  the  omni¬ 
directional  element  array  except  that  the  peak  response  at  =  0  *  and  9/  =  180  *  is  attenuated 
as  I  0,  I  increases.  This  is  due  to  the  X  axis  geophone’s  directional  response  and  because  the  Z 
axis  geophone  does  not  receive  the  signal  at  all.  At  0,  =  90  *  the  trace  is  clipped  at  -10  dB. 
Note  that  the  response  at  0,  =0  is  the  same  as  the  top  10  dB  of  the  bottom  panel  in  Figure  5.5. 

Figures  5.11  and  5.12  are  the  bearing  responses  of  the  omni-directional  and  geoplione 
element  arrays  with  the  source  at  9,  =80,  0,  =0  as  in  Figures  5.6.  The  response  along  the 
cone  swept  out  when  0,  0  *  and  0,  goes  from  0*  to  360*  is  hard  to  visualize.  However,  one 

can  see  that  the  (A,  =  0  *  traces  correspond  to  the  top  10  dB  of  the  plots  in  Figure  5.6.  Also, 
comparing  Figures  5.11  and  5.12  shows  that  the  source  direction  ambiguity  at  9/  =  120'  is 
reduced  by  the  geophone’s  directionality. 

The  array  bearing  response  is  even  harder  to  interpret  when  0,  ^  0  However,  a 
comparison  of  sparse  omni-directional  and  geophone  arrays  serves  to  emphasize  the  enhanced 
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Figure  5.9  Bearing  response  of  a  four  omniKiirccuonal  element  array  to  a  6.7  Hz  source  at  azimuth 
6,  =  0°  and  elevation  =  0^  The  array  clcmcnis  are  uniformly  spaced  at  100  meter  inter- 
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Figure  5.10  Bearing  response  of  a  four  element  array  of  orthogonally-oriented  gcophonc  triplets.  A  6.7 
Hz  source  is  located  at  azimuth  =  O*’  and  elevation  =  0°.  The  array  elements  are  uni¬ 
formly  spaced  at  100  meter  intervals. 
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Figure  S.ll  Bearing  response  of  a  four  omni-dircciionol  element  array  to  a  6.7  Hz  source  at  azimuth 
0#  =  and  elevation  =  0*.  The  array  elements  are  unifonniy  spaced  at  100  meter  inter¬ 
vals. 
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Figure  5.12  Bearing  response  of  a  four  element  array  of  orlhogonally-oricntcd  gcophonc  triplets.  A  6.7 
Hz  source  is  located  at  azimuth  0,  =  60®  and  elevation  =  0®.  The  array  elements  arc  uni¬ 
formly  spaced  at  100  meter  intervals. 
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capability  of  llie  gcoplioiie  array  to  resolve  source  direction.  Figure  .I  lS  shows  the  hearing 
response  of  the  omni-directional  element  array  with  the  source  at  0,  =  60,  <t>,  =  20.  The  array 
is  sparse  and  the  bearing  response  cojitains  numerous  grating  lobes  like  those  seen  in  Figure  5.7. 
The  line  array  directional  ambiguity  is  near  0i  =  120  '  as  long  as  the  look  elevation  is  small.  As 
the  look  elevation  angle  becomes  large  (positive  or  negative),  the  ambiguity  and  source  direction 
peaks  move  closer  together. 

Figure  5.14  shows  the  bearing  response  of  the  geophone  array  to  the  same  source.  The 
geophone  array  clearly  provides  an  enhanced  capability  to  determine  source  direction.  The 
level  in  the  direction  of  the  array  ambiguity,  di  =120  ’,  <i>,  =20,  is  5  dB  lower  than  in  the 
omni-directional  element  array.  The  strong  grating  lobe  seen  at  6,  «  205,  4),  «  —20  is  also 
.seen  in  Figure  5.8  and  5.13,  as  are  the  lower  level  grating  lobes.  There  has  been  significant  grat¬ 
ing  lobe  level  reduction,  although  there  remain  a  few  high  level  lobes  which,  in  the  presence  of 
noise,  might  be  confused  with  the  main  lobe. 


5.1.5.  Swallow  Float  Array 

As  a  last  simulation,  consider  the  array  formed  by  the  mid-water  column  floats  during 
the  September  1986  experiment.  The  element  spacing  is  about  1000  to  2200  meters,  so  the  array 
is  sparse,  and  the  elements  are  randomly  located  in  three  dimensions.  On  the  average,  sparse 
randomly-spaced  arrays  are  not  subject  to  grating  lobes,  but  may  have  side  lobes  whose  levels 
are  nearly  as  high  as  the  main  lobe.*®  Thus,  the  Swallow  float  array  bearing  response  may  be 
expected  to  contain  high  level  side  lobes. 

Figures  5.15  and  5.16  show  the  three-dimensional  bearing  response  of  a  four  element 
array  of  Swallow  floats  located  at  the  estimated  float  positions  for  record  1830  of  the  September 
1986  experiment.  The  two  Figures  contain  6.7  Hz  sources  at  azimuth  =  0  '  and  elevation 
4),  -  0  and  azimuth  0,  =60  '  and  elevation  4>,  =-10  ’,  respectively.  Both  Figures  contain 
numerous,  high  level  side  lobes,  many  only  1  dB  lower  than  the  main  lobe  level. 

The  bearing  response  for  small  elevation  angles  is  composed  of  ridges  which  are  narrow 
in  azimuth,  about  5 ' ,  and  broad  in  elevation,  about  20° .  This  vertical  ridge  structure  will  also 
be  seen  in  the  experiment  data. 

Other  prominent  features  are  a  strong  180°  ambiguity  and  a  weak  response  at  ±  90  ° 
from  the  source  direction.  Both  of  these  features  are  due  to  geophone’s  directionality,  and  they 
too  will  be  seen  in  the  experimental  data. 


5.2.  Application  to  Experiment  Data 

It  has  been  shown  that  the  proper  way  to  combine  orthogonally-oriented  geophone  sig¬ 
nals  is  to  multiply  by  the  look  direction  cosines  and  sum.  There  are,  however,  some  preliminary 
steps  which  must  be  taken  when  processing  experimental  data.  The  geophone  axes  are  not,  in 
general,  parallel  to  the  array  axes,  and  all  of  the  floats  do  not  have  the  same  time  base.  The 
steps  required  to  align  the  float  axes  with  the  array  axes  and  float  time  bases  with  each  other 
are  discussed  below,  after  which  the  experimental  data  are  beamformed.  The  experimental 
results  turn  out  to  look  much  like  noisy  versions  of  the  simulation  results. 


5.2.1.  Processing  the  Time  Series 

The  processing  steps  arc  listed  below  and  descril)ed  in  more  detail  in  subsequent  para¬ 
graphs. 


•  each  float’s  geophone  time  scries  is  operated  on  to  effectively  rotate  the  floats 

so  that  they  are  all  aligned  with  the  array  coordinate  system; 
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Figure  5.13  Bearing  response  of  a  four  omni-dircciional  clement  array  to  a  6.7  llz  source  at  azimuth 
0,  =  60°  and  elevation  =  20°.  The  array  elements  arc  uniformly  spaced  at  350  meter  inter¬ 
vals. 
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Figure  5.14  Bearing  response  of  a  four  element  array  of  onhogonally-orienicd  geophone  triplets.  A  6.7 
Hz  source  is  located  at  azimuth  0,  =  60*  and  elevation  6,  =  20*.  The  array  elements  are  uni¬ 
formly  spaced  at  350  meter  intervals. 
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Figure  5.15  Dcuring  response  of  a  four  clement  mray  of  Swallow  floats  located  at  estimated  float  posi* 
lions  for  record  1830  of  ihc  September  1986  experiment  Tlicre  is  a  6.7  Hz  source  at  azimuth 
0,  =  0®  and  elevation  =  0®.  Tlic  array  aperture  is  approximately  2000  meters. 
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Figure  5.16  Bearing  response  of  a  four  clement  array  of  Swallow  floats  located  at  estimated  float  posi¬ 
tions  for  record  1830  of  the  September  1986  experiment.  Tlicrc  is  a  6.7  Hz  source  at  azimuth 
0,  =  60®  and  elevation  =  -10®.  Tltc  array  aperture  is  approximately  2000  meters. 
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•  the  time  bases  of  all  but  one  of  the  floats’  time  series  are  shifted  to  align  them 
with  the  reference  float’s  time  base; 

•  each  axis  is  beamformed  separately;  the  X  axis  time  series  from  all  of  the 
floats  are  beamformed  together,  and  likewise  for  the  Y  and  Z  axis  time  series; 

•  the  beamformer  output  for  the  3  axes  are  multiplied  by  the  look  direction 
cosines  and  summed  to  form  a  time  series  for  each  look  direction; 

•  the  time  series  are  Fourier  transformed  and  a  calibration  curve  applied; 

•  the  levels  of  selected  spectral  peaks  are  selected  from  the  beam  spectra,  yield¬ 
ing  bearing  levels,  which  are  spectral  levels  as  a  function  of  bearing.  The 
bearing  level  is  analogous  to  the  magnitude  of  the  bearing  response  used  ear¬ 
lier. 


5. 2. 1.1.  Rotating  Float  Axes 

The  most  straightforward  approach  to  beamforming  Swallow  float  data  might  be  to 
combine  each  float’s  time  series  separately  to  form  a  beam  for  each  look  direction,  and  run  the 
beamformer  once  for  each  direction.  However,  this  approach  would  require  the  formation  of  a 
time  series  for  each  float  and  look  direction,  which  would  be  voluminous,  and  would  require  run¬ 
ning  the  beamformer  separately  for  each  look  direction,  which  is  inefficient  and  time  consuming. 

An  equivalent  approach  is  to  operate  on  each  float’s  time  series  to  align  their  axes  with 
the  array  axes,  beamform  each  axis  from  all  the  floats  separately,  and  combine  the  beamformer 
outputs  using  the  look  direction  cosines.  Using  this  latter  approach,  the  beamformer  is  only  run 
once  for  each  axis  and  only  one  time  series  per  look  direction  is  produced.  This  is  the  method 
used  in  this  thesis. 

Rotating  the  floats  to  align  their  axes  with  the  array  axes  can  be  understood  by  refer¬ 
ring  to  Figure  5.17,  which  is  a  plan  view  of  a  freely-drifting  float  within  the  coordinate  system 
established  for  localizing  the  floats.  The  direction  of  the  float’s  y  axis  is  denoted  Ojoj,,  and  is 
shown  in  the  Figure  to  be  about  120’  magnetic.  The  float’s  z  axis  is  vertical  and  parallel  to  the 
Z  axis  of  the  coordinate  system. 

The  time  series  from  the  float  geophones  which  are  parallel  to  the  float’s  x  and  y  axes, 
are  denoted  x{n)  and  y(n).  They  must  be  combined  to  obtain  new  gcophone  time  scries  in  the 
directions  of  the  array  X  and  Y  axes,  call  them  x'(n}  and  y'(n).  By  inspection, 

i'(n)  =  |y(n)  cos  (280  "  -  Oflo.t)  +*(«)  sin  (280  ’  -Ofloat)] 

and 

y'(n)  =  l^y(n)  sin  (280  *  -  -j(n)  cos  (280  *  -Osoat)]  (5-20) 

Equations  (5.20)  show  how  each  float’s  geophone  time  series  can  be  combined  to  effectively 
rotate  the  float  axes  to  align  them  with  the  array  axes.  The  angle  is  the  float  heading, 

which  the  float  measures  at  the  start  of  each  record.  The  directions  of  the  coordinate  system 
axes  are  determined  by  plotting  the  positions  of  the  bottomed  floats  on  a  nautical  chart. 


5. 2. 1.2.  Aligning  Float  Time  Bases 

The  floats  have  different  time  bases  because  their  internal  clocks  are  not  identical,  i.e., 
the  clock  rates  vary  slightly  from  float  to  float.  Consider  floats  i  and  /.  Denote  time  according 
to  float  j’s  clock  as  T^.  According  to  float  j,  record  n-)-l  starts  at 
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Figure  6.17  Plan  view  of  a  freely-drifting  Swallow  float,  whose  heading  is  denoted 
a  fifth  ^he  coordinate  system  defined  by  the  positions 
of  bottomed  floats  0  and  1 


Tj=45n  (5.21) 

which  is  45  n  of  float  y’s  seconds.  In  (21]  it  is  shown  that  once  the  float  temperatures  stabilize, 
the  relationship  between  elapsed  time  according  to  a  float,  and  true  elapsed  time  since  syn¬ 
chronization,  is  linear.  Let  true  elapsed  time  be  t.  Then 


Ty(0  =  a, 


(5.22) 


where  jy  and  ay  are  constants  unique  to  each  float.  Note  that  jy  is  dimensionless  and  the  units 
of  Oy  are  time. 

Equation  (5.22)  can  be  solved  for  true  time  /  in  terms  of  Ty,  jy  and  ay,  and  if  gj  and  Oy 
are  known,  the  float  time  base  can  be  aligned  to  true  time.  However,  it  is  shown  in  [21]  that 
the  py’s  and  Oy’s  cannot  be  measured  directly.  It  is  further  shown  in  [21]  that  reciprocal  path 
travel  times  can  be  combined  to  estimate  the  clock  rate  and  offset  differences,  (j,-  —  jy)  and 
(a,-  —  ay).  Therefore,  the  approach  is  to  pick  one  of  the  floats  as  the  reference  time  base  and 
align  all  of  the  other  floats  to  it. 

Solving  (5.22)  for  true  time  t  and  using  (5.21)  yields  the  true  elapsed  time  at  the  start 
of  record  n-f-1  aboard  float  j. 
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*») 


(5.23) 


Changing  only  the  subscript  yields  the  true  elapsed  time  at  the  start  of  record  n+1  aboard  float 

I . 


<(i.  n)  =- 


Ti  —  a,-  45  n  —  a,- 


9i 


9i 


(5.24) 


The  difference  between  (5.24)  and  (5.23)  gives  the  difference  between  the  time  bases  of  the  two 
floats  in  true  time. 

^  «(».  y>  n)  =  t{j,  n)  -  t{i,  n) 


45  n 


9) 


45  n  —  a,- 
9i 


=  45  n 


9i 


i 

—  (true  time  dif.  between  time  b&ses) 

9) 


(5.25) 


Figure  5.18  depicts  the  time  bases  of  floats  i  and  j  and  that  of  true  time  when  float  I’s 
clock  rate  is  somewhat  faster  than  float  j's.  The  horizontal  axes  in  the  top  panel  are  time 
according  to  float  i,  time  according  to  float  and  true  time.  All  are  increasing  to  the  right, 
which  means  that  l(y,  n)  ><(i,  n)  and  A  <(i,  j,  n)  >0. 

The  event  indicated  by  l(e)  in  Figure  5.18  occurs  simultaneously  at  the  two  floats. 
However  the  time  of  the  event  according  to  float  i  is  later  than  the  time  of  the  event  according 
to  float  j. 


[T..(c)-T,(n)]>[r,.(e)-T»] 

If  the  time  series  are  simply  aligned  using  the  record  start  times,  as  shown  in  the  bottom  panel 
of  Figure  5.18,  then  the  event  does  not  occur  simultaneously  in  the  two  time  series  but  occurs 
first  at  float  j .  In  order  for  the  two  events  to  add  coherently,  they  must  occur  simultaneously  in 
the  time  series  presented  to  the  beamformer.  One  or  the  other  of  the  float  time  bases  must  be 
shifted  to  align  it  with  the  other. 

Let  float  j  be  chosen  as  the  reference  time  base,  so  that  float  I’s  time  base  is  to  be 
shifted  A  l(i,  j ,  nl  seconds.  To  find  the  amount  of  the  shift  in  terms  of  float  I’s  clock,  denoted 
A  Ti  I A  t{i,  i,  n)  ,  use  equation  (5.22)  and  substitute  A  t(i,  j,  n)  for  (: 

A  Ti  I A  <(«,  j,  rj)j  =  gi  |a  <(i,  ;,  n)j  +  a. 


Using  (5.25), 


’  ' 

' 

[a  <(«,  y,  «)]  =  ?, 

45  n 

J _ 1_ 

9i  9,  , 

9i  9i 
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event 


Figure  5.18  Comparison  of  the  lime  bases  of  floats  i  and  j  and  mie 

time.  The  upper  panel  shows  their  relationship  when  the  floats 
were  deployed.  The  bottom  panel  shows  the  effect  of  simply 
aligning  the  record  start  limes. 
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=  45  n 


9i  -  9} 
Oj 


+  a, - ay 

9} 


(5.26) 


Note  that  if  jy  =  j,  and  oy  =a,-,  which  is  to  say  that  the  two  time  bases  are  the  same,  then 
A  T,  =0  as  expected. 

It  is  shown  in  [2l]  that  the  values  of  are  very  close  to  1,  such  that  if  we  write 
=  1  +  €,  then  £  is  of  order  10”*.  Thus  (5.26)  can  be  simplified  using  the  Taylor  series  expan¬ 
sion 


!{x  4-  a)  =/(i)  -t-/'(x)  -t-  •  ■  • 


of 


Using 


/'(x)=-L  and  r(x)=-T 

the  expansion  is 

1 _ 1  _  g  ^ 

X  +  a  X  x®  X® 


Substituting  x  =  1  and  a  =  £,  to  first  order  in  (, 


_1_ 

9i 


1 

1  -t-<, 


1  -£, 


Therefore,  to  first  order  in  £, 


9i  -  9j  ^  1  -b  £,■  -  1  -  £y 

9i  l+f; 


(«.• 


-£y)  (1  -£y)  «  £,•  -  £y  =  -  Jy 


(5.27a) 


and 


ii.  =  J-lii 

9i  1  +  «y 


(1  -b  £,)  (I  —  £y)  »  1  -b  £<  —  fy  1  -b  9i  “  Jy 


(5.27b) 


Putting  equations  (5.27)  into  (5.26)  yields 

A  T,  «  45  n  [y,  -  ffy]  -b  a,-  -  [l  -b  -  ?y]  oy  =  ^45  n  -  gy  j  -  jy  j  -b  o.-  -  gy  (5.28) 
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It  is  shown  in  [21]  that 


I  ay  1  <1  second  <C45  n  seconds,  n  >1 

Equation  (5.28)  therefore  becomes 

A  T,  «  45  n  I?,-  -  ffy  j  -I-  [fl,-  -  ay  j 


(5.29) 


(5.30) 


It  is  now  shown  that  A  l{i,  j,  n)  changes  negligibly  over  the  course  of  one  record  so 
that  it  is  sufficient  to  realign  the  time  series  once  for  each  record.  To  see  this,  use  (5.25)  to  find 
the  change  in  A  <(■)  over  the  course  one  record. 


A  <(«,  j,  n+l)  -  A  t{i,  j,  n) 


J_ 

Qi 


To  first  order  in  £, 


1  -  £y  -  1  -)-  £,  =  £,  -  £y 


SO  that 


A  t(»,  j,  n-l-1)  —  A  t(i,  j,  n)  «  45  £,•  —  «y  45  •  10“*  seconds  «  0.5  msec. 


(5.31) 


Equation  (5.30)  thus  shows  how  to  calculate  how  much  float  i’s  time  base  must  be 
shifted  at  the  start  of  each  record  in  order  to  align  it  with  the  corresponding  record  in  float  j’s 
geophone  data.  It  requires  knowing  the  differences  between  the  clock  rates  and  offsets  in  the 

two  floats,  —  jy  ■  and  |a,  —  ay  |.  Reference  [2l)  shows  how  the  differences  between  the  clock 


rates  and  offsets  are  calculated,  and  contains  values  for  them  for  all  float  pairs  during  the  Sep¬ 
tember  1986  deployment.  The  procedure  for  aligning  the  float  time  bases  is: 


•  chose  one  float  as  a  reference  whose  time  base  will  not  be  shifted; 

•  using  measured  values  for  (j,-  —  ?y)  and  (a,-  —  ay),  and  equation  (5.30),  calcu¬ 
late  how  much  each  record  of  the  other  float’s  geophone  data  must  be  shifted, 
and  in  what  direction; 

•  interpolate  and  shift  the  three  geophone  time  series  for  all  but  the  reference 
float. 


6. 2. 1.3.  Remaining  Processing  Steps 

Having  aligned  the  gcophone  time  series  in  direction  and  time,  they  can  now  be 
delayed  and  summed  to  listen  selectively  in  particularly  directions.  A  time  delay  beamformer  is 
applied  to  the  time  series  from  each  axis  of  all  the  floats,  producing  three  time  series  for  each 
look  direction.  Let  6,(1,  0t,  4>i]  denote  the  time  series  for  look  direction  (^(,  4>i),  obtained  by 
combining  the  geophones  parallel  to  the  i  axis.  The  three  sets  of  time  series  are  then  multiplied 
by  the  look  direction  cosines  and  summed. 


61 


R.  L.  Culver 


*(t.  <t>i)  =  **(<.  ^1.  <^1)  cos(^,)  cos(^,)  +  b^{t,  01,  <t>,)  sin(^,)  cos(^,) 


+  b,{t,  <l>t)  sin(^,) 


(5.32) 


The  time  series  for  each  look  direction  is  Fourier  transformed  using  the  same  length 
transform  (40.96  seconds)  and  window  used  for  the  spectra  in  Appendix  B. 


B{f,0„4>,)^F 


(5.33) 


After  applying  a  calibration  curve,  the  levels  of  selected  spectral  lines  may  be  picked  from  all  of 
the  beam  spectra  to  obtain  the  spectral  level  as  a  function  of  bearing  or  bearing  level.  The 
bearing  level  is  analogous  to  the  magnitude  of  the  bearing  response  which  was  used  in  Section 
5.1.  It  indicates  the  directionality  of  sound  at  particular  frequencies. 


5.2.2.  Horizontal  Bearing  Levels 

The  initial  effort  is  to  beamform  only  in  the  horizontal  direction,  that  is,  for  between 
0°  and  360”  and  (t>i  =0”.  This  is  considered  a  reasonable  approach  because  the  bearing 
responses  in  Figures  5.15  and  5.16  show  that  the  horizontal  and  low  elevation  angle  bearing 
responses  look  much  alike  as  long  as  the  source  arrival  angle  is  small,  and  according  to  Figure 
3.10,  the  vertical  arrival  angles  are  small  for  many  records  at  the  frequencies  of  interest. 

Horizontal  bearing  levels  for  a  20  record  (15  minute)  period  beginning  with  record  1800 
are  shown  in  Figure  5.19.  The  upper  5  dB  of  the  bearing  levels  at  9,  10.3,  12.2  and  13.5  Hz  are 
shown.  The  azimuth  angle  is  degrees  relative  to  magnetic  north  and  the  look  elevation  is  0” . 

The  four  plots  in  Figure  5.19  contain  numerous,  high  level  side  lobes  and  a  180*  ambi¬ 
guity  much  like  Figures  5.15  and  5.16.  The  azimuth  angle  of  the  highest  lobe  is  not  constant  for 
more  than  a  few  records  (3  to  5  minutes)  at  a  time.  In  general,  however,  the  source  bearing 
appears  to  be  near  either  090”  or  270’  magnetic.  A  means  of  enhancing  the  main  lobe  and 
reducing  the  level  of  the  side  lobes  is  needed  in  order  to  better  resolve  the  source  direction. 


5.2.3.  Enhancing  Bearing  Levels 

Figure  5.20  contains  the  result  of  averaging  the  bearing  levels  of  five  consecutive 
records  at  each  frequency.  The  entire  110  record  period  selected  for  analysis  is  shown  in  these 
plots.  Averaging  appears  to  smooth  the  bearing  levels,  so  that  the  source  direction  is  approxi¬ 
mately  the  center  of  a  rounded  mass  of  peaks.  Gross  changes  in  the  approximate  source  bearing 
during  the  110  record  period  are  approximately  the  same  at  all  frequencies.  For  example,  the 
source  azimuth  increases  about  20’  between  records  1780  and  1800.  This  indicates  that  all  of 
the  signals  originate  from  the  same  source. 

As  another  approach  to  enhancing  the  bearing  levels.  Figure  5.21  contains  the  product 
of  the  9,  10.3  and  12.2  Hz  bearing  levels.  If  different  frequency  signals  arrive  at  the  array  from 
the  same  direction,  their  main  lobes  will  line  up  but  their  side  lobes,  whose  spacing  is  frequency 
dependent,  will  not.  Figure  5.21  shows  that  multiplying  bearing  levels  together  thins  the  peaks 
out  considerably,  but  the  180’  ambiguity  persists.  This  implies  that  the  signal  arrival  angles 
are,  to  some  extent,  frequency  dependent. 

Figure  5.22  summarizes  the  horizontal  direction  beamforming  results.  The  top  panels 
show  the  highest  bearing  level  during  the  110  record  period.  Comparing  the  top  panel  of  Figure 
5.22a  with  Figure  3.10  shows  that  the  gross  rise  in  signal  level  evident  in  the  individual  float 
spectral  lines  is  also  seen  in  the  bearing  levels.  This  is  true  for  the  other  frequencies  as  well. 
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Figure  S.19a  Spectral  level  versus  azimuthal  bearing  at  9  liz  from  four  freely-drifting  Swallow  floats  dur¬ 
ing  records  1800  -  1819  of  the  September  1986  experiment  (IS  minutes  duration).  The  look 
elevation  angle  is  zero  (Itorizontal). 
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Figure  5.19b  Spectral  level  versus  azimuthal  bearing  at  10.3  Hz  from  four  freely-drifting  Swallow  floats 
during  records  1800  -  1819  of  the  September  1986  experiment  (IS  minutes  duration).  The 
look  elevation  angle  is  zero  (horizontal). 
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Figure  5.19c  Spectral  level  versus  azimuthal  bearing  at  12.2  Hz  from  Tour  freely-drifting  Swallow  floats 
during  records  1800  -  1819  of  the  September  1986  experiment  (15  minutes  duration).  The 
look  elevation  angle  is  zero  (horizontal). 
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Figure  5.20a  Averaged  spectral  levels  versus  a/imuihal  bearing  at  9  Hz  from  four  freely-drifting  Swallow 
floats  during  records  1758  -  1868  of  the  September  1986  experiment  (82.5  minutes  duration). 
The  look  elevation  angle  is  zero  (horizontal).  Each  trace  is  the  average  of  five  records. 
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Figure  5.20b  Averaged  speciral  levels  versus  azimuthal  bearing  at  10.3  llz  from  four  ftccly-drifling  Swal 
low  floats  during  records  1758  -  1868  of  the  September  1986  experiment  (82.S  minutes  dura 
tion).  The  look  elevation  angle  is  zero  (horizontal).  Each  trace  is  the  average  of  five  records. 
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Figure  5.20c  Averaged  spectral  levels  versus  ar-imulhal  bearing  at  12.2  Hz  from  four  freely-drifting  Swal¬ 
low  floats  during  records  1758  -  1868  of  the  September  1986  experiment  (82.5  minutes 
duration).  The  look  elevation  angle  is  zero  (horizontal).  Each  trace  is  the  average  of  five 
records. 
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Figure  5.20d  Averaged  spectral  levels  versus  azimuthal  bearing  at  13.5  Hz  from  four  frccly-drifling  Swal¬ 
low  floats  during  records  1758  -  1868  of  the  September  1986  experiment  (82.5  minutes 
duration).  The  look  elevation  angle  is  zero  (horizontal).  Each  trace  is  the  average  of  five 
records. 
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Figure  5.21  Product  of  9,  10.3  and  12.2  Hz  spccual  levels  versus  azimuthal  bearing  from  four  freely- 
drifting  Swallow  floats  during  records  1800  -  1819  of  the  September  1986  experiment  (IS 
minutes  duration).  The  look  elevation  angle  is  zero  (horizontal). 
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The  middle  panels  in  Figure  5.22  show  the  exact  frequency  at  which  the  peaks  occur. 
Comparing  them  with  Figure  3.9  shows  that  changes  in  signal  frequency  which  are  apparent  in 
the  individual  float  spectra  also  appear  in  the  bearing  levels.  The  bottom  panels  in  Figure  5.22 
show  the  apparent  direction  to  the  source  based  upon  the  azimuth  of  the  highest  bearing  level. 
The  source  direction  jumps  between  90“  and  approximately  270° ,  indicating  a  180"  ambiguity. 


5.2.4.  Three-Dimensional  Bearing  Levels 

Appendix  D  contains  bearing  levels  for  record  1802  for  azimuth  angles  between  150“ 
and  240“  relative  to  the  array  X  axis  and  elevation  angles  between  -20°  and  20“  .  The  arrival 
angles  during  this  record  were  seen  in  Figure  3.11  to  be  nearly  horizontal.  The  three- 
dimensional  bearing  levels  display  the  same  vertical  ridge  structure  seen  in  Figures  5.15  and 
5.16.  The  bearing  levels  at  different  frequencies  all  display  main  lobes  at  azimuth  angle 
6i  «  187  but  the  elevation  angles  vary  between  0“  and  15  “ . 

Appendix  D  also  contains  bearing  levels  for  record  1844  for  azimuth  angles  between 
150“  and  240“  (relative  to  the  array  X  axis)  and  elevation  angles  between  -60“  and  60  “.  This 
record  was  selected  because  the  magnitude  of  the  arrival  angle  at  9  Hz  is  seen  in  Figure  3.11  to 
be  greater  that  30“  for  three  of  the  floats.  The  9  Hz  bearing  level  indicates  that  considerable 
energy  is  arriving  at  large  positive  and  negative  elevation  angles  during  record  1844.  This  indi¬ 
cates  that  the  signal  is  propagating  to  the  array  over  multiple  paths.  Energy  at  10.3  Hz  is  con¬ 
centrated  between  -15“  and  15“  elevation.  The  other  two  signals,  at  12.2  and  13.5  Hz,  are 
arriving  at  -45“  to  20“  elevation  and  -35”  to  40“  elevation,  respectively.  We  see  that  the 
vertical  angle  arrival  varies  considerably  with  frequency. 


5.2.5.  Beamformer  Gain 

One  method  of  evaluating  how  well  signals  are  being  co-phased  is  to  compare  the 
highest  bearing  level  with  the  levels  at  the  individual  floats.  If  the  beamformer  is  perfectly  co¬ 
phasing  each  float’s  signals,  then  the  highest  bearing  level  will  be  equal  to  the  coherent  sum  of 

f  -2  ' 

the  signals.  The  level  in  dB  of  a  sinusoidal  signal  with  amplitude  x,  is  10  log 


coherent  sum  of  N  equal  strength  signals  having  the  same  phase  is 


The 


10  log 


{Nx?  _ 


20  log  N  +10  log 


(5.30) 


which  shows  a  gain  of  20  log  N  over  the  level  of  one  signal.  When  there  are  4  elements,  the  gain 
is  12  dB. 

When  the  signals  do  not  all  have  the  same  amplitude,  as  is  the  case  with  the  Swallow 
floats,  but  tlhey  are  co-phased,  the  coherent  sum  is 

10  log  (xi -f- 12  +  *3  +  2^4)^  (coherent  sum)  (5-31) 


where  z,  is  the  signal  amplitude  at  float  i.  When  the  signals  have  approximately  the  same 
amplitude,  the  coherent  sum  is  something  less  than  12  dB  higher  than  the  level  of  the  individual 
signals. 

Figure  5.23  compares  the  highest  horizontal  bearing  levels  to  the  coherent  sum  of  the 
spectral  levels  of  each  float’s  x  and  y  axis  geophones.  The  coherent  sum  has  been  subtracted 
from  the  highest  bearing  level.  Also  shown  is  the  incoherent  sum  of  the  x  and  y  axis  geophone 
levels.  If  all  of  the  signals  are  of  equal  strength,  the  incoherent  sum  is 
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Figure  5.22a  Hor^ntal  beamforming  results.  DirccUon.  frequency  and  level  of  the  highest  horizontal 
bearing  level  for  the  9  Hz  spectral  line. 
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Figure  5.22c  Horizontal  bcamforming  results.  Direction,  frequency  and  level  of  the  highest  horizontal 
bearing  level  for  the  12.2  Hz  spectral  line. 
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10  log  (  N  —  )  =  10  log  N  10  log 


o 


so  that  the  signal  gain  with  N  =  4  is  10  log  4  =  6  clB.  The  coherent  sum  lias  also  been  sub¬ 
tracted  from  the  incoherent  sum  in  Figure  5.23,  and  we  see  that  the  incolierent  sum  is  roughly  6 
dB  lower  than  the  coherent  sum. 

The  highest  bearing  levels  are,  on  the  average,  about  3  dB  lower  than  the  coherent 
sum  and  about  3  dB  above  the  incoherent  sum.  This  indicates  that  the  bearing  level  gain  is 
about  9  dB.  The  signals  are  apparently  being  co-phased  well  enough  for  the  highest  bearing 
level  to  be  higher  than  the  incoherent  sum,  but  are  not  being  perfectly  co-phased,  in  that  the 
highest  bearing  level  is  lower  than  the  coherent  sum. 

The  three  x’s  on  the  first  three  plots  in  Figure  5.23  mark  the  highest  three-dimensional 
bearing  levels,  normalized  by  the  coherent  sum  of  the  horizontal  and  vertical  geophone  spectral 
levels.  Including  the  vertical  channel  does  not  change  the  conclusion  regarding  bearing  level 
gain. 


5.2.6.  Correlating  Geophone  Signals 

In  a  somewhat  separate  effort  to  determine  the  direction  to  the  source,  geophone  sig¬ 
nals  from  pairs  of  floats  were  correlated.  The  correlation  of  two  broadband  signals  is  a  spike  at 
a  particular  time  delay.  Narrowband  signals  have  a  sinusoidal  correlation  function.  The  exper¬ 
imental  data  contain  ship-generated  signals,  which  are  approximately  sinusoidal,  but  it  was 
hoped  that  they  had  sufficient  bandwidth  for  the  correlation  functions  to  yield  unique  time 
delays.  The  delays  for  all  float  pairs  might  then  be  used  to  estimate  the  source  bearing  from 
each  pairs  of  floats. 

The  time  series  from  the  4  floats  were  Fourier  transformed  and  pairwise  multiplied 
together  to  obtain  their  cro.ss  power  spectra.  The  cross  power  spectra  were  then  divided  by  the 
power  spectrum  of  one  of  the  pair  to  obtain  the  transfer  functions  of  the  channels  separating 
each  float  pair.  The  channel  transfer  functions  were  inverse  Fourier  transformed  to  obtain  the 
channel  impulse  responses.*®  In  order  to  restrict  the  correlation  to  frequencies  originating  from 
the  source,  only  the  17  bins  centered  at  each  frequency  were  used. 

The  correlation  results  are  shown  in  Appendix  E.  They  contain  many  peaks,  indicating 
that  there  was  insufficient  bandwidth  to  the  source  lines  for  the  correlation  method  to  yield  use¬ 
ful  information  about  source  direction. 
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Figure  S.23a  Highest  horizontal  bearing  level  and  incoherent  sum  of  the 
9  Hz  (upper  panel)  and  13.5  Hz  Gower  panel)  spectral  levels 
of  each  float’s  x  and  y  axis  geophones.  Both  are  normalized 
by  the  coherent  sum  of  all  floats’  horiTontal  spectral  levels. 

X’s  maiit  the  highest  normalized  3-d  bearing  levels. 
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Figure  5.23b  Highest  horizontal  bearing  level  and  incoherent  sum  of  the 

10.3  Hz  (upper  panel)  and  12.2  Hz  (lower  panel)  spectral  levels 
of  each  float’s  x  and  y  axis  geophones.  Both  arc  normalized 
by  the  coherent  sum  of  all  floats’  horizontal  spectral  levels. 
X’s  mark  the  highest  normalized  3-d  bearing  levels. 
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6.  Analysis  and  Conclusions 

The  beamforming  results  shown  in  Chapter  5  indicate  that  different  frequency  signals 
emanating  from  the  same  source  arrive  at  the  array  at  different  azimuth  and  elevation  angles, 
and  at  different  strengths.  Also,  because  the  floats  are  so  widely  dispersed  and  the  array  aper¬ 
ture  so  large,  the  azimuth  and  elevation  arrival  angles  at  a  particular  frequency  vary  consider¬ 
ably  from  float  to  float.  This  indicates  that  signals  are  propagating  to  the  array  by  frequency- 
dependent  paths,  and  are  not  propagating  across  the  array  as  plane  waves.  Since  the  beam- 
former  calculates  delays  to  co-phase  signals  propagating  across  the  array  as  plane  waves,  this 
may  explain  why  the  peak  bearing  levels  are  not  a.s  high  as  the  coherent  sum. 

The  effect  of  improperly  co-phasing  the  time  series  (due  to  non-plane  wave  propaga¬ 
tion)  can  be  simulated  by  feeding  the  beamformer  erroneous  float  positions,  causing  it  to  calcu¬ 
late  an  erroneous  set  of  time  delays.  The  effect  of  different  strength  signals  can  be  simulated  by 
introducing  shading  errors.  The  difference  in  float  spectral  levels  at  the  frequencies  of  interest 
was  found  to  be  14  dB  or  less,  except  when  an  anomaly  such  as  clipping  occurs.  Figure  6.1 
shows  the  effect  of  shading  errors  alone  on  the  horizontal  bearing  response  of  the  four  Swallow 
float  array  during  record  1830  with  a  10.3  Hz  sovjrce  at  azimuth  202°  ,  elevation  0°  .  Shading 
errors  uniformly  distributed  between  0  and  14  dB  have  been  introduced  into  each  trace.  It  is 
apparent  that  this  type  of  error  by  itself  does  not  produce  results  like  those  seen  in  the  experi¬ 
mental  data,  e.g.  Figure  5.19. 

Figure  6.2  shows  the  bearing  response  of  the  same  array  to  the  same  source,  but  here 
position  errors  have  been  combined  with  shading  errors.  The  position  errors  used  in  this  simula¬ 
tion  are  mean  zero  and  have  a  standard  deviation  of  20  meters.  Shading  errors  are  uniformly 
distributed  between  0  and  14  dB.  This  Figure  looks  much  like  the  experimental  bearing  level 
plots  in  Figure  5.19. 

The  experimental  bearing  levels  in  Figure  5.19  indicate  that  the  source  is  either  to  the 
east  or  to  the  west  of  the  array.  The  rough  underwater  terrain  to  the  east  makes  it  more  likely 
that  the  source  is  actually  east  of  the  array,  so  that  the  high  bearing  levels  at  azimuth  270° 
magnetic  are  due  to  the  180°  ambiguity.  Relatively  steep  underwater  cliffs  begin  about  5  miles 
east  of  the  array  and  rise  to  quite  shallow  water,  approximately  100  meters  deep,  about  14 
miles  to  the  east.  These  rough  features  may  cause  sound  to  be  scattered  in  many  directions,  so 
that  from  the  array,  the  source  direction  appears  to  change  erratically  with  time. 

The  shallow  water  east  of  the  array  also  offers  a  plausible  explanation  for  the  absence 
of  the  4.5  Hz  fundamental  frequency  in  the  Swallow  float  data.  The  cutoff  frequency  may  be 
higher  than  4.5  Hz  over  some  portion  of  the  transmission  path  between  the  source  and  the 
array.  To  calculate  the  water  depth  required  to  block  transmission  of  a  4.5  Hz  signal,  assume  a 
3  layer  model  consisting  of  a  constant  sound  speed  water  layer  bounded  from  above  by  an  air 
half-space  and  from  below  by  a  sediment  half-space.  As  shown  in  Chapter  9  of  Medwin  and 
Clay,^®  the  equation  which  sets  the  minimum  depth  for  a  particular  frequency  to  propagate  is 


(m  -  l)ff  +  <t>,  +  <t>t  < 


2  ff  A 


1  - 


12 


1/2 


(6.1) 


where  m  is  the  mode  number,  2  is  the  phase  angle  of  the  reflection  coefficient  at  the  upper 
(air-water)  boundary,  2  is  the  phase  angle  of  the  reflection  coefficient  at  the  lower  (water- 
sediment)  boundary,  h  is  the  water  depth,  cq  and  c,  are  the  sound  speeds  in  the  water  and  sedi¬ 
ment,  respectively,  and  X  is  the  wavelength  of  sound  at  the  frequency  of  interest.  Assuming  an 
ideal  pressure  release  air-water  interface,  2  =  jt.  Assuming  an  ideal  rigid  bottom,  =0.  Set 

m  equal  to  1  because  the  first  mode  must  propagate  for  any  mode  to  propagate.  The 
wavelength  at  4.5  Hz,  assuming  a  sound  speed  of  1500  meters/sec,  is  approximately  330  meters. 
Substituting  these  values  into  equation  (6.1)  yields 
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Bearing  Response,  Record  1830,  10.3  Hz.  Source  @  202 
No  Position  Error,  <  14  dB  Shodfng  Error 


<S>  T3 


:  Azimuth  Angle,  Degrees  CCw  from  X  Axis 

Figure  6.1  Horizontal  bearing  response  of  the  four  Swallow  float  array  during  record  1830  to  a  10.3  Hz 
source  at  azimuth  202®,  elevation  0®.  showing  the  effect  of  0  to  14  dB.  uniformly  distributed 
shading  errors. 


Power  Densi'ty,  dB 
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Figure  6.2  Horizontal  bearing  response  of  iltc  four  Swallow  float  array  during  record  1830  to  a  10.3  Hz 
source  at  azimuth  202°,  elevation  0°,  showing  the  effect  of  shading  errors  and  mean  zero, 
normally  distributed  position  ''rrors  (standard  deviation  =  20  meters). 
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and  the  condition  on  h  is 

£o 

«i 

If  Cl  Cq,  then 

A  >  75  meters  c,  Cq 

If  Cl  =  1.2  Cq,  then 

A  >  135  meters  Cj  =  1.2  Cq 

Thus  depending  upon  the  values  chosen  for  Ci  and  Cq,  the  water  depth  would  need  to  be  shal¬ 
lower  than  about  150  meters  in  order  to  block  a  4.5  Hz  signal.  We  saw  in  Figure  3.2  that  this 
condition  is  met  to  the  east  of  the  array.  It  is  not  necessary  that  the  source  be  in  shallow 
water,  only  that  the  transmission  path  include  shallow  water. 

The  variation  in  vertical  arrival  angle  seen  in  Figure  3.11  may  have  something  to  do 
with  the  radiation  pattern  of  the  source.  If  the  9  Hz  line  set  are  blade  pass  harmonics,  they  are 
radiated  by  the  tip  of  each  blade  of  the  ship’s  propeller  as  it  passes  the  hull.  The  sound  radi¬ 
ates  strongly  outward  and  downward.  If  the  other  lines  at  10.3  and  12.2  Hz  are  auxiliary 
machinery  tonals,  they  tend  to  radiate  primarily  outward  rather  than  downward.  Thus  we 
might  expect  the  9,  13.5  Hz,  etc.  lines  to  arrive  at  higher  angles  than  the  10.3  and  12.2  Hz  sig¬ 
nals. 
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APPENDIX  B 


Geophone  Power  Spectra 


Float  5,  X  Axis  Geophone  Power  Spectra 
Records  1810  -  1821  and  the  12  Record  fiverogi 
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Comparison  of  Kalman  and  Least  Squares  Filters 
for  Locating  Autonomous  Very  Low  Frequency  Acoustic  Sensors 


R.  L.  Culver  and  W.  S.  Hodgkiss 


Marine  Physical  Laboratory 
Scripps  Institution  of  Oceanography 
San  Diego,  CA  92152 


1.  Introduction 

Under  Office  of  Naval  Research  sponsorship,  the  Marine  Physical  Laboratory  has 
designed  and  built  self-contained,  freely-drifting  Swallow  floats  which  can  measure  ambient 
ocean  noise  in  the  1  to  20  Hz  band  for  up  to  25  hour  periods.  The  floats  are  ballasted  for  neu¬ 
tral  buoyancy  at  a  specified  depth.  The  deployment  of  several  floats  forms  a  freely-drifting 
array  of  noise  sensors  whose  outputs  may  be  combined  after  the  experiment  with  a  beamformer. 
[1] 

Float  locations  must  be  known  in  order  to  beamform  their  outputs.  The  floats  gen¬ 
erate  and  receive  acoustic  pulses  and  thereby  measure  float  to  surface  and  float  to  float  travel 
times.  A  serie.s  of  signal  processing  routines  have  been  developed  to  estimate  float  position  from 
these  travel  time  measurements.  This  paper  compares  two  localization  methods: 

•  A  generalized  least  squares  filter  estimates  current  position  using  the  current  set  of 
travel  time  measurements.  The  filter  uses  only  that  set  of  measurements  and  an 
estimate  of  the  mean  square  measurement  error. 

•  A  Kalman  filter  estimates  current  position  by  weighting  the  current  set  of  meas¬ 
urements  against  the  previous  position  estimate  propagated  forward  in  time.  Posi¬ 
tion  estimates  are  propagated  forward  in  time  using  equations  of  motion  developed 
for  the  floats.  The  filter  also  uses  estimates  of  the  measurement  error  variance 
and  float  accelerations. 

Computer  simulation  is  used  to  compare  filter  performance  under  several  deployment 
scenarios.  The  effects  of  increasing  the  randomness  of  float  movement  and  the  measurement 
error  are  also  investigated.  Results  show  that  the  Kalman  filter  performs  better  than  the  least 
squares  filter  as  long  as  the  floats  are  only  subjected  to  small  magnitude  accelerations  between 
measurements.  This  is  true  even  when  there  are  large  measurement  errors.  Neither  filter  was 
found  to  be  sensitive  to  relatively  major  changes  in  deployment  geometry  as  long  as  the  sound 
speed  profile  is  known  exactly.  However,  deploying  the  floats  in  a  vertical  array  does  degrade 
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the  performance  of  both  filters. 

MPL  has  been  deploying  Swallow  floats  since  1982.  A  number  of  articles  and  technical 
reports  have  been  written  describing  various  aspects  of  the  project.  [2-11)  The  reader  is  referred 
to  articles  cited  in  the  reference  list  for  details  about  Swallow  float  construction  and  functioning 
not  found  in  this  paper. 


2.  Travel  Time  Measurements  and  Mean  Square  Measurement  Error 

In  this  section  we  discuss  the  travel  time  measurements  made  by  the  floats  and  show 
data  from  a  September  1986  Swallow  float  experiment.  We  describe  how  the  measurements  are 
processed  to  remove  errors  caused  by  differences  among  float  clock  rates  and  estimate  mean 
square  measurement  error.  The  resulting  float  to  float  and  float  to  surface  travel  time  and 
mean  square  error  estimates  are  the  inputs  to  the  Kalman  and  least  squares  filters. 

The  general  Swallow  float  hardware  configuration  is  shown  in  Figure  1.  Suspended 
about  1  meter  below  the  glass  sphere  is  an  acoustic  transducer  which  generates  and  receives  8 
kHz  tone  bursts.  Pulses  10  milliseconds  long  are  transmitted  by  the  floats  in  a  pre-programmed 
sequence.  A  different  float  transmits  every  45  seconds.  When  there  are  N  floats  deployed,  each 
float  transmits  every  45- ;V  seconds.  Twelve  floats  are  normally  deployed,  and  each  float  nor¬ 
mally  transmits  every  9  minutes. 

The  floats  are  listening  whenever  they  are  not  transmitting.  They  receive  pulses 
transmitted  by  other  floats  as  well  as  reflections  of  their  own  pulses.  They  ordinarily  detect 
reflections  from  the  ocean  surface,  and  often  from  the  ocean  bottom. 

The  float  pulse  detection  circuit  outputs  an  analog  signal  indicating  whether  or  not  a 
pulse  is  detected.  The  signal  is  sampled  by  the  float  at  1  kHz  {1000  times  per  second)  and  a  bit 
is  set  to  1  if  a  pulse  is  present,  0  otherwise.  Groups  of  eight  detection  circuit  samples  are 
inspected  by  the  float  for  the  presence  of  range  pulses  (bits  which  are  set).  Bytes  which  contain 
bits  set  are  stored  internally  along  with  their  time  according  to  the  float’s  internal  clock. 

After  a  deployment,  the  stored  pulse  arrival  information  may  be  used  to  reconstruct 
the  sampled  pulse  detection  circuit  output.  We  can  display  this  information  directly  by  con¬ 
verting  it  to  a  raster  image.  Figure  2  is  one  such  plot,  taken  from  the  September  1986  deploy¬ 
ment  of  12  Swallow  floats  in  the  Pacific  ocean  50  miles  west  of  San  Diego,  Ca.  The  vertical 
scale  of  the  plot  has  been  changed  from  travel  time  (in  seconds)  to  range  (in  meters)  using  the 
approximate  speed  of  sound,  1500  meters/second.  The  plot  shows  detections  made  by  float  4 
during  records  when  float  3  transmitted.  The  depth  of  both  floats  was  about  1000  meters  and 
water  depth  was  2000  meters.  [8] 

The  line  of  very  short  dashes  running  across  the  upper  part  of  the  plot  mark  the 
^  arrival  of  the  direct  path  pulse  from  float  3  at  float  4.  These  dashes  are  about  15  meters  (10 
msec)  in  length  and  are  separated  by  12  records.  The  much  longer  vertical  lines  running  across 
the  plot  at  greater  range  mark  the  arrival  of  pulses  which  bounced  off  the  ocean  surface  and 
then  traveled  to  the  float.  These  lines  are  150  to  600  meters  (100  to  400  msec)  long  and  are  also 
separated  by  12  records.  The  arrival  of  surface  bounce  pulses  is  spread  out  in  time  because  of 
scattering  by  the  sea  surface. 

A  computer  program  has  been  developed  to  detect  the  upper  edge  of  each  dash  or  line, 
which  corresponds  to  when  the  pulse  first  arrives.  The  program  is  referred  to  as  the  edge  detec¬ 
tor  [5]  in  that  it  outputs  the  time  corresponding  to  each  pulse’s  leading  edge.  Applying  the  edge 
detector  to  the  pulse  arrival  information  shown  in  Figure  2,  we  obtain  the  time  series  shown  in 
Figure  3.  Both  Figures  are  plotted  on  the  same  scale  so  that  they  may  be  compared  to  verify 
edge  detector  performance. 

Noise  in  the  edge  detector  output  is  evidenced  by  the  bad  point  at  about  record  1850. 
Smaller  magnitude  measurement  noise,  too  small  to  be  visible  on  this  scale,  is  also  present  in  the 
edge  detector  output.  [12] 
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The  pulse  arrival  time  according  to  the  edge  detector  minus  the  time  when  the  pulse 
was  transmitted  is  the  travel  time  measurement.  Each  float  determines  the  pulse  transmission 
time  according  to  its  own  internal  clock.  However,  the  clock  rates  vary  slightly  among  the 
floats  so  that  reciprocal  path  travel  time  measurements  must  be  interpolated  and  averaged  to 
obtain  the  true  interfloat  travel  time.  There  are,  in  fact,  several  processing  steps  which  the 
measurement  time  series  must  pass  through  before  being  fed  into  the  Kalman  or  least  squares 
filters.  These  steps  are  shown  in  Figure  4.  The  dotted  line  shows  how  simulated  travel  time 
measurements  are  introduced  in  order  to  validate  the  processing  and  filtering  under  different 
conditions. 

A  one-step  forward  adaptive  linear  predictor  |5]  is  used  to  reject  noisy  points  (like  the 
one  evident  in  Figure  3)  and  fill  gaps  which  sometimes  occur  in  the  travel  time  series.  The 
linear  predictor  must  be  told  the  maximum  allowable  difference  between  the  predicted  value  and 
a  data  point  for  the  data  to  be  considered  valid.  Predictor  output  is  monitored  to  ensure  that 
it  is  functioning  properly  and  if  necessary,  the  maximum  difference  parameter  is  adjusted. 

We  must  interpolate  the  travel  time  measurements  before  averaging  reciprocal  path 
measurements  because  the  measurements  are  not  made  simultaneously.  They  are  separated  in 
time  by  at  least  45  seconds,  and  by  as  much  as  N/2  *  45  seconds,  where  N  is  the  number  of 
floats  deployed.  A  cubic  spline  interpolation  is  used. 

Reciprocal  path  interfloat  travel  times  are  then  averaged.  The  averaged  interfloat 
travel  times  and  float  to  surface  travel  times  are  ready  for  input  to  the  least  squares  or  Kalman 
filter. 

An  estimate  of  the  mean  square  measurement  error  is  obtained  by  taking  the  mean 
square  of  the  difference  between  reciprocal  path  travel  time  measurements  minus  a  least  squares 
linear  fit  to  the  difference.  [12]  The  mean  square  measurement  error  estimate  is  used  by  both 
filters  to  weight  the  measurements,  that  is,  give  greater  weight  to  those  with  smaller  mean 
square  error. 


3.  The  Least  Squares  Filter 

Our  problem  is  to  use  the  travel  time  measurements  and  measurement  error  estimates 
just  described  to  estimate  float  positions.  The  least  squares  filter  operates  on  one  set  of  travel 
time  measurements  and  produces  a  set  of  float  positions  which  minimizes  the  square  of  the 
differences  between  the  measurements  and  the  interfloat  travel  times  calculated  using  the  posi¬ 
tion  estimate.  A  generalized  least  squares  estimator  is  employed,  allowing  greater  weight  to  be 
ipplied  '■"rtain  measurements.  The  meas^irement  error  estimates  determine  the  weights,  with 
greater  weight  applied  to  measurements  whose  error  is  thought  to  be  lower. 

The  derivation  of  the  least  squares  estimator  and  the  loss  function  linearization 
method  presented  here  closely  follows  that  given  by  Sorenson.  [13]  We  start  by  identifying  the 
n‘*  measurement  vector 


=C' 


'!  2 


'I  3 


^N-{  N 


(3.1) 
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where  f,"  is  the  averaged  travel  time  estimate  between  floats  i  and  j ,  f"  is  the  round  trip  travel 
time  estimate  between  float  i  and  the  surface,  and  C  is  a  matrix  whose  diagonal  elements  are 
the  estimated  average  sound  speeds  for  the  paths  over  which  the  travel  times  were  measured. 


‘‘nh  n 


The  last  N  elements  of  (3.2),  which  arc  divided  by  2,  correspond  to  round  trip  travel  times. 
The  measurement  vector  is  related  to  the  float  position  vector  according  to 

Z„^h{X,)  +  V, 


A',  is  the  vector  of  float  positions  at  time  n 


Vi  ••  •  Vn 


and  x"  is  the  x  coordinate  of  float  1  at  time  n,  etc.  The  origin  of  the  coordinate  system  is  at 
the  ocean  surface  directly  above  float  1,  and  z  is  positive  downward.  The  axes  are  rotated  so 
that  float  2  is  in  the  xy  plane. 

/i(X„)  indicates  a  nonlinear  transformation  on  A',  from  position  space  to  measurement 
space  (ranges  and  depths). 


j 

[(xj  -xr)2+(y.?  -yr)2+(^5  -Xf)']' 

J 

[(^.^  +(y;  -yrf +(^; 


ma;)  = 


J 

[(*N  ~  +  (jw  ~  1/N-I  f  +  ~  ^,V-1  )^| 


102 


SIO  Reference  88-16 


The  first  group  of  elements  of  A  (AT,)  are  the  interlloat  ranges  and  the  second  group  are  the  float 
to  surface  distances,  or  depths. 

V.  n  equation  (3.3)  is  the  measurement  error  or  noise  vector 


Cl  2*^1  2 


^N-l  N^N-i  N  C,V, 


<^N^N 


(3.6) 


whose  elements  are,  or  are  made  to  be,  mean  zero.  The  subscripts  of  the  v’s  are  the  same  as 
those  of  the  e  ’s  and  r ’s  described  below  (3.1).  The  covariance  of  V,  is  R ,  where 


E 


(3.7) 


VVe  assume  that  the  noise  is  uncorrelated  from  measurement  to  measurement,  an  assumption 
which  has  been  demonstrated  to  be  reasonable.  (12)  R  is  thus  diagonal,  with  elements  (12) 


(3.8) 


where  f,"  and  e'y  are  identified  in  (3.1)  and  (3.2),  and  and  <7^  are  the  variances  of  the  uncer¬ 
tainties  in  e,y  and  r,y,  respectively. 

The  least  squares  estimate  of  AT,,  given  measurements  Z„,  is  defined  as  the  value 
which  minimizes  the  following  loss  function: 


/ 


LS  = 


^ZjR,  Z,  -2ZlR,  h(X,)+h(Xy  /?,  k{X,) 


(3.9) 


Necessary  and  sufficient  conditions  for  X,  to  minimize  lis  are  that  dli^fdX  =0  and 
is  positive-definite.  Thus  to  determine  the  least  squares  estimate  X,  form 


dX 


=  -  2Zl  R. 


dhjX) 

dX 


+  2h{Xy  R, 


dhjX) 

dX 


=  2 


h{X^]-Z. 


R. 


dh{X) 

dX 


=  0 


(3.10) 


Equation  (3.10)  says  that  the  residual  h(X„)—Z^  must  be  orthogonal  to  the  columns  of 
dh{X)/dX.  If  h{X)  is  a  linear  function  of  X 

A(XJ  =  //  X. 


and 


dh[X,] 

dX 


=  // 
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then  (3  10)  represents  a  linear  equation  in  X  with  solution 


Xls  = 


R  H 


//^  R  Z 


h{X) 


linear  in  x 


(3.11) 


However,  when  h(X)  is  a  nonlinear  function  of  X,  as  in  (3.5),  (3.10)  represents  a  system  of  non¬ 
linear  algebraic  equations  for  which,  in  general,  there  is  no  closed  form  solution.  As  a  result,  we 
must  linearize  the  loss  function  and  solve  the  problem  numerically.  A  numerical  procedure  will 
not  in  general  provide  an  exact  solution  to  (3.10).  Our  hope  is  that  it  provides  an  improved 
estimate  of  X,  so  that  repeated  application  may  yield  a  satisfactory  estimate. 


Three  numerical  methods  were  considered  for  estimating  X„  given  Z,  and  R„.  In  what 
follows,  we  identify  k  as  a  step-size  parameter  which  mu.st  be  chosen.  To  simplify  the  notation, 
define 


A  ^(X) 
'  ^  OX 


X  =  x 


(3  12) 


In  order  of  increasing  complexity,  the  methods  considered  are  the  steepest  descent  procedure 


-A(X,)] 


(3.13) 


=  A'.  -1-  k,  HT  R 

the  Gauss  procedure 

X.^  =i,  k,  [hT  R  //,]■*  Hi  R  [z  -  M^.)]  (3.14) 

and  the  Newton  procedure 

X,^  =  X,  +  k,  j/y.^  R  -  S  [^y  -  h,(X,)]  %  ^(^.)|  '  Hi  R  [z  -  A(X,)]  (3.15) 


Here  the  subscript  i  is  the  iteration  index. 

.Newton’s  method  was  chosen  for  our  use  because  it  provides  the  most  rapid  conver¬ 
gence  in  the  vicinity  of  the  minimum  value  of  lig.  (13)  Newton’s  method  does  entail  more  com¬ 
putation  than  the  other  methods,  but  that  is  not  a  problem  for  us  at  this  time. 

We  now  derive  the  Newton  procedure  to  see  what  assumptions  are  made.  We  start 
with  an  initial  estimate  of  X„,  which  we  call  X,,  a  vector  of  interfloat  range  and  float  depth 
measurements  Z„,  and  a  covariance  matrix  R„.  We  want  to  find  a  new  estimate  of  X„  which 
reduces  the  value  of  the  loss  function  identified  in  equation  (3.9).  First  we  use  a  Taylor  s.eries  to 
expand  the  loss  function  in  the  vicinity  of  X„,  keeping  terms  of  order  2  or  less  in  (X,  —  X„). 


/(XJ  «  /(XJ 


iMl 

dX 


X.  -  X , 


(X,  -  XJ-t-^(X, -X.)’ 


d^nx) 

dX^ 

(X.  -X.08.16) 


Necessary  and  sufficient  conditions  for  the  loss  function  to  be  minimized  are  that  d^l /dX^  is 
positive  definite  and  that  the  gradient  of  the  loss  function  vanish.  The  gradient  of  (3.16)  is 


dX 


-‘(X) 


dX 


a^/(X) 

5X^ 


X.  -X. 


{X,  -X.) 


(3.17) 
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Defining 


M2Q 

dX 


(3.18) 


dV) 

dX^ 


d 

T 

X.  -A'.  ^  dX 

dX 

=  •*  ’n 


iv. 


(3.19) 


and 


(3.20) 


the  requirement  that  (3.17)  vanish  may  be  written  compactly: 

9n+KP.=  0  (3.21) 

g„  and  P„  are  vectors;  iV„  is  a  matrix.  and  are  calculated  using  Equation  (3.21)  is 

solved  using  the  least  squares  method  discussed  below  to  obtain  P„.  We  then  substitute  it  into 
(3.20)  to  obtain  a  new  estimate  of 

(old)X„-f-P,  t  (new)X„  (3.22) 

We  must  now  check  our  new  value  of  X  ^  against  some  criterion  and  decide  if  it  is  a 
satisfactory  estimate  of  X„.  The  criterion  used  is  discussed  below.  If  X ^  is  satisfactory,  we  are 
finished  and  can  go  on  to  find  an  estimate  for  the  float  positions  at  the  time  of  the 

(n  +  l)‘*  measurement.  If  it  is  not,  we  recalculate  g„  and  using  the  new  value  of  X^,  obtain 
a  new  least  squares  solution  for  P„,  and  obtain  still  a  newer  value  of  X„  using  (3.22).  This  pro¬ 
cedure  may  be  repeated  until  a  satisfactory  estimate  of  X„  is  found. 

In  summary,  the  estimate  of  is  found  from  Z„  and  P,  using  the  following  steps; 


Summary  of  Least  Squares  Procedure 

Step  1:  Pick  X„. 

Step  2:  Determine  if  X„  is  a  satisfactory  estimate  of 
X„.  If  so,  terminate  the  search,  increment  n 
and  go  to  Step  1.  If  not,  go  to  Step  3. 

Step  3;  Calculate  and  using  equations  (3.18)  and 
(3.19)  and  obtain  a  least  squares  solution  for 
Pn- 

Step  4:  Calculate  a  new  value  for  X ^  using  equation 
(3.22)  and  return  to  Step  2 
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Fixing  the  Coordinate  System 

We  fix  the  coordinate  system  so  that  the  origin  lies  at  the  surface  directly  above  float 
1,  and  float  2  lies  in  the  xz  plane.  Tljis  is  accomplished  by  specifying  that  the  coordinates  of 
floats  1  and  2  in  the  initial  estimate  are  such  that  float  1  is  on  the  z  axis  and  float  2  is  in 

the  xz  plane. 


Then  P-  is  constrained  such  that 


=  V?  =  !/?  =  0 


1  0  0 
0  1  0 
0  0  0 
0  0  1  1  fo 
0  0  01  =  |0 
0 


0  0  0 


(3.23) 


(3.24) 


Since  each  new  estimate  is  the  sum  of  the  old  X  „  and  P,,  the  constraint  (3.24)  ensures  that 
equation  (3.23)  continues  to  be  satisfied. 


The  Least  Squares  Solution  for  P„ 

Obtaining  the  least  squares  solution  to  (3.21)  subject  to  the  constraint  (3.24)  is  accom¬ 
plished  using  a  method  recommended  by  Lawson  and  Hanson.  (14)  We  first  derive  an  (N  —  l)  ■  3 
dimensional  unconstrained  least  squares  problem,  solve  the  derived  problem,  and  finally 
transform  the  solution  of  the  derived  problem  to  obtain  a  solution  to  the  original  constrained 
problem.  The  QR  method,  which  involves  repeated  application  of  the  Householder  transforma¬ 
tion,  is  used  to  solve  the  derived  least  squares  problem. 


Search  Termination  Criterion 


The  search  termination  criterion  is  used  to  decide  if  A"  „  is  a  satisfactory  estimate  of 
A„.  We  specify  that  the  search  will  stop  when  the  mean  magnitude  of  the  residual  changes  by 
less  than  1  meter  for  10  consecutive  values  of  A„.  The  mean  magnitude  of  the  residual  is 


mean  magnitude  of  the  residual  = 


I 

M 


(3.25) 


where  M  is  the  number  of  elements  in  the  vector  Z.  When  10  consecutive  values  of  (3.25) 
change  by  less  than  1  meter,  we  stop  the  search  and  accept  the  current  value  of  X,. 


10 


16 


SIO  Reference  88-16 


4.  The  Kalman  Filter 

The  Kalman  filter  represents  an  attempt  to  obtain  a  better  estimate  of  the  float  posi¬ 
tions  by  taking  advantage  of  additional  information  which  may  be  available.  This  additional 
information  consists  of. 

•  a  model  of  how  the  floats  propagate  through  the  water, 

•  an  estimate  of  the  magnitude  of  the  accelerations  to  which  the  floats  are  subjected 
while  they  are  deployed,  and 

•  an  estimate  of  the  variance  of  the  error  in  the  initial  estimate  of  the  float  positions. 

When  these  additional  pieces  of  information  are  incorporated  in  the  Kalman  filter,  it  can  per¬ 
form  better  than  the  least  squares  filter  under  conditions  which  we  will  determine. 

The  form  of  the  Kalman  filter  is 


Here  X,  is  the  current  estimate  of  the  float  positions  and  velocities.  is  the  previous  esti¬ 

mate.  The  vectors  and  A(<I>  X  ,_i)  are  unchanged  from  the  definitions  in  (3.1)  and  (3.5), 
respectively.  The  vector  X„  is  now  twice  as  long  as  it  was  in  the  least  squares  filter  derivation 
because  it  now  includes  float  velocities  as  well  as  positions. 

<1>  is  called  the  transilion  matrix  and  it  propagates  the  float  positions  forward  in  time. 
K  is  called  the  Kalman  gain  matrix.  Its  magnitude  determines  how  much  the  current  float  posi¬ 
tion  estimate  is  influenced  by  the  current  measurements  as  opposed  to  relying  primarily  upon 
the  previous  estimate  propagated  forward  in  time. 

There  are  many  books  and  articles  dealing  with  the  Kalman  filter.  [15-20]  The  deriva¬ 
tion  given  here  closely  follows  Sorenson.  [17] 

We  start  with  a  model  of  the  floats  as  they  propagate  through  the  ocean.  For  simpli¬ 
city,  consider  for  the  moment  only  one  float.  Its  position  at  time  n  is  (z",  y",  z")  and  its  velo¬ 
city  components  are  (i*,  y",  i").  The  position  coordinates  and  velocity  components  are  com¬ 

bined  into  the  vector  X,,  called  the  state  vector. 

xn  [*"  y-  y"  ^*]  (4.2) 

We  model  the  float’s  movement  through  the  water  with  the  following  linear  difference  equations: 

x"+'  =  I-  +  i"  T  +  a;^  z*-”  =  i*  -I- o;t 

=y"  +y-r y"-^  =  y"  -b  a;7 


=  z*  -b  i"  r -t- af—  =  i"  -b  a;r  (4.3) 


107 


R.  L.  Culver 


where  T  is  the  elapsed  time  between  n  and  n  +  I.  The  accelerations 


(i.4) 


are  unknown,  but  have  zero  mean 


E{yl)  =0 


and  known  covariance 


E 


A,  /\l]  a  QA 


(4.5a) 


(4.5b) 


Equation  (4.5b)  is  part  of  the  new  information  which  is  used  by  tlie  Kalman  .liter  but  not  by  the 
least  squares  filter.  VVe  define  the  transition  matrix 


<1>  ^ 


1  0  0  T  0  0 
0  10  0  TO 
0  0  1  0  0  T 
0  0  0  1  0  0 
0  0  0  0  1  0 
0  0  0  0  0  1 


(4.6) 


and 


H 

2 


0 


r  t 


0 


0  0 


T  0 
0  T 
0  0 


0 

0 

'Jl 

2 

0 

0 

T 


Using  (4.6)  and  (4.7),  we  can  combine  equations  (4.3)  into  one  vector  equation: 


(4.7) 


=<l>  X,  +r  A,  (4.8) 

Equation  (4.8)  is  referred  to  as  the  dynamic  or  plant  equation  or  process.  We  can  expand  X,,  4>, 
and  r  so  that  equation  (4.8)  describes  the  motion  of  all  of  the  floats.  Equation  (4.8)  is  a  new 
equation  which  was  not  used  by  the  least  squares  filter  and  which  will  be  incorporated  into  the 
Kalman  filter. 

The  exact  value  of  the  state  at  time  0  is  unknown.  However,  we  assume  that  its 
expected  value 
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E(Ao)— Xq  (4.9a) 

and  covariance  matrix 

Po  ^  E  [(Xo  -  Xo)  (Xo  -  i'o)'' ]  (4.9b) 

are  known.  Equation  (4.9b)  is  also  new  information  which  is  used  by  the  Kalman  filter. 

The  Kalman  filter  makes  use  of  the  plant  equation  (4  8)  and  the  measurement  equation 
identified  in  (3.3).  Our  goal  is  to  determine  an  estimate  of  A',  which  is  a  linear  combination  of 
the  estimate  and  the  measurements  Z,  and  which  minimizes  the  expected  value  of  the  sum 
of  the  squares  of  the  errors  in  the  estimate.  That  is, 

E  [(X,  -X,)^  (X,  -X,)|  is  minimized  (4.10) 

The  condition  (4.10)  is  equivalent  to  minimizing  the  trace  of  P,,  the  covariance  of  the  estimate 
error. 

E  [(X.  -X.)^  (X„  -Xj]  =  trace  |  E  [(X.  -X.)  (X.  -X.)^]  =  trace  P.  (4.11) 


We  know  that  the  state  of  the  floats,  that  is  their  positions  and  velocities,  evolves 
according  to  (4.8).  Because  the  accelerations  A  are  mean  zero,  one  wpuld  not  expect  them  to 
influence  the  state  from  one  measurement  to  the  next.  Thus  given  X  it  is  reasonable  to 

predict,  in  the  absence  of  any  measurement,  that 

X,  =<t>X,_i  (before  any  measurement)  (4  12) 

The  measurement  at  time  n  can  be  used, to  modify  this  estimate.  Based  upon  (4.12)^one  would 
expect  the  n‘*  measurement  to  be  A(4>  X,_,).  The  difference  between  Z,  and  A(<i>  X,_,)  is  the 
new  information,  sometimes  called  the  innovation,  contained  in  the  measurement.  We  define  the 
estimate  of  X„  to  be 

X,  =<i>x,_,-(-a:.[z, -A(<I>X,^)]  (4.13) 

where  the  matrix  X,  will  be  determined  to  satisfy  the  condition  (4.10). 

The  error  in  X,  is,  using  (4.13)  and  (3.3), 

X„  ^X,  -X. 

=  <1>X._^+X,[z.  -/,(^x,_,)]-x. 

«  <!>  X,_,  +  X,  [z.  -  M«*>  X,^}]  -  <t>  x._, 
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We  need  to  linearize  li(X)  in  order  to  simplify  (4.14).  We  can  use  a  Taylor  series 
h(X)  about  its  nominal  value,  which  we  call  ,V. 


Using  (4.15)  to  simplify  h{^  A’«-i)  -  A(<1>  X,_,),  we  get 


(X-X) 


h{X,)-h{X,) 


y  n)  +  gy 


X  -  .V* 


{X,  -X,) 


h{X\)  + 


dhlX) 

dX 


h{X\)-h(X<)  + 


^h{X) 

ax 


x-x. 


dhjX) 

ax 


x.i,  ix,  -X'J 


It  seems  reasonable  in  light  of  (4.12)  to  choose  for  our  nominal  values 

.V.  =A". 

This  reduces  (4.16)  to 


h{X,)-h(X,) 


dh(X) 

ax 


Putting  this  result  into  (4.14)  yields 


A,  «  <i>X^^-K. 


ah(x) 

ax 


+x,  V, 


-K, 


dh(X) 

ax 


X  -♦  X 


t>  +  /c,  K. 


I  -K. 


ah{x) 

ax 


X.*X^_, 


t>X,^+I<,  V, 


Now  we  form  the  covariance  of  the  estimate  error 


E 


>.  -V/]  k 


no 


(4.14) 
to  expand 

(4.15) 

{X„  -X\) 

(4.16) 

(4.17) 


(4.18) 

(4.19) 
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Putting  (4.18)  into  (4.19),  using  (3.7),  and  assuming  that  X  and  V  are  uncorrelated,  we  can 
eventually  obtain  [17] 


P.-\l-K  (4.20) 


where  we  have  identified 


K  ^  +  T  A  r’ 


To  simplify  the  notation,  define 


A  mx) 

*=  dX 


(4.22) 


Then  (4.20)  can  be  written 


P.  =  M,  -  K,  //.  M„  -  M,  +  K,  [//.  A/.  +  P.  ]  (4.23) 

The  term  A/„  +  P*  j  is  symmetric  and  non-negative  definite,  so  it  can  be  written  as  the 

product  of  a  matrix  5,  and  its  transpose  (i.e.  a  matrix  square  root). 


S.  M„  H^  +  R 


(4.24) 


Substituting  (4.24)  into  (4.23)  and  postulating  the  existence  of  an  array  A  in  order  to  complete 
the  square,  we  obtain 

P„  =  M,  -  K,  H„  M„  -  M,  Hi  Kl  +  Sl  Kl 


=  M„+  K,S„  Si  Kl  -  K„  S^  Al  -  A,  5.  AT.  +  A,  A, 


-  K„  //,  A/„  -  A/,  Hi  Kl  +  K,  5„  Al  +  A,  5„  K„  -  A.  Al 


=  A/„  4-  [k„  5,  -  A„ )  (P„  5,  -  A, )%  K„  [s,  Al  -H„M,]+  [a„  Sl  -  Af,  HI\kI  - 


>1,  Al 


Now  define  A,  so  that 


S,  Al  -  H„  A/„  =  0 


A,  =M„  HHs;^f 


(4.25) 
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When  5,  is  singular,  the  pseudo  inverse  is  used  instead  of  Using  (4.25),  the  error  covari¬ 
ance  becomes 

P.  =  A/.  +  [k,  5.  -  )  (/C.  5.  -  A,  Al  (4.26) 

Only  the  product  term  of  (4.26)  involves  /v, .  The  product  of  a  matrix  and  its  transpose  is 
positive-definite.  Thus  the  trace  of  P  is  minimized  by  choosing 

(4.27) 


The  Kalman  gain  matrix  is  then 

K,  =  M„  Hi  (5,'')^  5,-* 


=  A/,  Hi  (Si  5J-‘ 

=  A/.  Hi  [//.  A/.  Hi +  /i„ 

In  terms  of  the  covariance  of  the  estimate  error  is 


P,  P,  Af. 


(4.28) 


(4.29) 


The  matrix  inversion  lemma  [17]  may  be  used  to  rewrite  K„  in  terms  of  the  covariance  matrix 


K,  =  P,  Hi  P,-* 


(4.30) 


where 


P„=M,-  Hi  H„  M,  Hi  +  R,\  //,  M, 


(4.31) 


Equations  (4.13),  (4.21),  (.4.30),  and  (4.31)  comprise  the  Kalman  filter  estimator.  In 
summary,  we  obtain  the  estimate  X„  using  the  following  steps: 


Fixing  the  Coordinate  System 

We  want  the  coordinate  system  to  be  fixed  as  in  the  least  squares  filter.  To  do  this,  we 

specify 


•  X Q  such  that  float  1  is  on  the  z  axis  and  has  zero  horizontal  velocity,  and  float  2  is  on 

the  xy  plane  and  has  zero  velocity  in  the  y  direction; 
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Summary  of  Kalman  Filter  Procedure 

Step  1: 

Pick  Xq  and  Pq.  Set  n  =  1.  Go  to  step  4. 

Step  2: 

Using  (4.21),  calculate  the  estimate  error  co- 
variance  before  the  measurement. 

Step  3: 

Using  (4.31),  calculate  the  estimate  error  co- 
variance  after  the  measurement. 

Step  4: 

Calculate  /i(4>  X,_,)  using  (3. .5)  and  //„  using 
(4.22). 

Step  5: 

Calculate  the  Kalman  gain  using  (4.30). 

Step  6: 

Calculate  the  position  and  velocity  estimates 
using  (4.13). 

Step  7: 

Increment  n  and  go  to  Step  2. 

•  Pq  such  that  the  errors  in  the  initial  estimates  of  float  I’s  *  and  y  coordinate  and  velo¬ 
city  and  the  errors  in  the  initial  estimates  of  float  2’s  y  coordinate  and  velocity  are 
known  to  be  zero. 

•  Qq  such  that  float  1  is  known  to  experience  zero  horizontal  acceleration  and  float  2  is 
known  to  experience  zero  acceleration  in  the  y  direction. 


5.  Float  Location  Simulations 

We  now  proceed  to  test  the  functioning  of  the  travel  time  measurement  processing  dis¬ 
cussed  in  Section  2  and  programs  implementing  the  least  squares  and  Kalman  filters  derived  in 
Sections  3  and  4.  Test  cases  are  developed  by  simulating  float  motion  on  the  computer  and  cal¬ 
culating  the  distance  between  the  floats  at  periodic  intervals  to  simulate  travel  time  measure¬ 
ments.  The  measurements  are  then  processed  and  fed  to  the  filters.  The  resulting  position  esti¬ 
mates  are  compared  to  the  actual  float  positions.  We  calculate  the  rms  position  error  for  each 
float  and  average  over  all  of  the  floats  to  obtain  a  performance  measure. 

Figure  5  shows  the  steps  taken  to  simulate  travel  time  measurements.  We  first  estab¬ 
lish  the  initial  float  positions.  Four  different  deployment  geometries  are  simulated  in  order  to 
investigate  their  effect  on  filter  performance.  We  then  simulate  float  motion,  propagating  them 
through  space  by  applying  mean-zero  random  accelerations  in  the  horizontal  directions  and 
oscillatory  motion  in  the  vertical  direction,  simulating  the  influence  of  internal  waves  in  the 
ocean.  Our  experience  deploying  and  tracking  Swallow  floats  in  the  ocean  has  given  us  a  feel 
for  how  much  the  floats  can  be  expected  to  move  horizontally  when  they  are  drifting  freely  and 
the  kind  of  vertical  motion  amplitudes  and  frequencies  which  can  be  expected  to  occur.  We 
increase  the  horizontal  and  vertical  excitation  (called  proctss  noise)  in  two  simulations  in  order 
to  see  its  effect  on  filter  performance. 
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Simulated  travel  time  measurements  are  derived  from  llie  float  positions  by  calculating 
the  distance  between  each  pair  of  floats  and  dividing  by  the  speed  of  sound.  Of  course  we  must 
introduce  errors  into  the  measurements  in  order  to  make  the  simulation  realistic.  We  add  errors 
(usuall  called  measurement  noise)  to  these  measurements  and  delete  a  few  points  to  simulate 
data  dropouts  because  they  occasionally  occur  in  real  travel  time  measurements.  The  measure¬ 
ment  noise  variance  in  two  simulations  wa.s  varied  in  order  to  determine  its  effect  on  filter  per¬ 
formance.  No  errors  were  introduced  in  the  sound  speed  estimates. 

The  four  float  deployment  geometries  used  in  the  simulations  are  depicted  in  Figure  6. 
Cases  1,2,6  and  7  used  our  normal  Swallow  float  deployment  configuration  in  which  we  deploy 
3  floats  to  the  ocean  bottom  (2000  meters  in  this  case)  and  5  (or  more)  floats  to  somewhere  in 
the  water  column  (1000  meters  depth  in  ihi.s  case),  'fhe  bottomed  floats  form  a  triangular  base 
for  locating  the  freely  drifting  floats.  The  mid-water  column  floats  form  an  approximately 
planar  array  which  is  used  to  measure  ambient  ocean  noise 

Case  1  establishes  the  baseline  for  evaluating  how  the  filters  respond  to  increased  noise 
or  special  deployment  geometries  Both  the  process  and  measurement  noise  are  low  in  Case  1. 
Cases  2,  6,  and  7  test  the  filters’  performance  when  the  process  and  measurement  noise  are 
increased.  The  proces.s  noise  is  raised  by  a  factor  of  100  in  Crises  6  and  7.  The  measurement 
noise  is  raised  by  a  factor  of  2  in  Cases  2  :ind  6. 

Case  3  tests  the  filters’  performance  when  the  mid-water  column  floats  drift  out  of  the 
net  formed  by  the  bottomed  floats  The  bottomed  floats  are  the  same  as  in  Case  1,  but  the 
planar  array  has  been  moved  considerably  outside  the  triangle  formed  by  the  bottomed  floats. 
Low  process  and  measurement  noise  are  use<l  for  this  simulation. 

Case  4  tests  the  effect  of  one  of  the  bottomed  floats  failing.  The  remaining  bottomed 
floats  and  planar  array  are  the  same  as  in  Case  I,  and  low  process  and  measurement  noise  are 
used. 

Case  5  tests  our  localization  capability  when  the  mid-water  column  floats  are  deployed 
in  a  vertical  line  array  rather  than  a  planar  array.  This  geometry  was  used  during  the  Sep¬ 
tember  1987  experiment  conducted  in  the  Pacific  Ocean  approximately  260  miles  west  of  Point 
Conception,  California.  In  this  simulation,  the  three  bottomed  floats  were  put  at  about  4000 
meters  depth  and  5  other  floats  were  spaced  between  700  and  2200  meters  depth.  Low  process 
and  measurement  noise  were  used  for  this  simulation. 


Table  5.1.  Float  Location  Sinmulations. 


Case 

Process  Noise 

Measurement  Noise 

T?sts  the  Effect  of: 

1 

low 

low 

filter  performance  baseline 

2 

low 

high 

increased  measurement  noise 

3 

low 

low 

large  array  drift 

4 

low 

low 

losing  a  bottomed  float 

5 

low 

low 

vertical  line  array 

6 

high 

high 

increased  process  noise 

7 

high 

low 

increased  process  and  measurement  noise. 

The  least  squares  filter  requires  an  initial  estimate  of  the  float  positions  and  the  Kal¬ 
man  filter  requires  an  initial  estimate  of  the  float  positions  and  velocities.  We  used  estimates 
which  were  nearly  identical  to  the  actual  float  positions  and  velocities  for  our  simulations.  This 
was  done  in  order  to  isolate  their  response  to  changes  in  process  noise,  measurement  noise,  or 
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deployment  geometry  from  problems  associated  with  locking  on  to  the  float  trajectories.  Train¬ 
ing  the  Kalman  filter  is  a  separate  problem  not  addressed  here. 

The  results  of  the  simulations  are  summarized  in  Table  5.2  and  plotted  in  Figures  7 
and  8.  Figure  7  compares  the  results  of  Cases  1,  2,  6  and  7  and  shows  the  effects  of  increased 
process  and  measurement  noise.  Figure  8  compares  the  results  of  Cases  1,  3,  4  and  5  and  shows 
the  effect  of  different  float  deployment  geometries. 


Table  5.2.  Summary  of  Float  Location  Simulation  Results. 

Average  RMS  Error,  meters 

Case 

Kalman  Filter 

Least  Squares  Filter 

0.9 

1.4 

1.6 

2.9 

0  7 

1.4 

■■ 

0.8 

1.5 

5 

1.2 

2.3 

6 

1.7 

1.7 

7 

2.3 

2.7 

Conclusions 

The  simulation  results  lead  one  to  the  following  conclusions; 

•  When  the  process  and  measurement  noise  are  low,  (Case  1),  the  error  in  the  least 
squares  position  estimate  is  about  50%  greater  than  that  in  the  Kalman  filter  estimate. 

•  Increasing  the  process  noise  increases  the  error  in  the  Kalman  filter  estimate  but  has 
no  effect  on  the  least  squares  estimate.  As  the  process  noise  increases,  the  Kalman 
filter  estimate  error  approaches  the  least  squares  estimate  error.  At  sufficiently  high 
process  noise,  the  two  filters  perform  comparably. 

•  When  the  process  noise  is  low,  increasing  measurement  noise  has  a  stronger  negative 
effect  on  the  least  squares  estimate  error  than  on  the  Kalman  filter  estimate  error. 
This  becomes  less  true  as  the  process  noise  is  increased. 

•  Moving  the  mid-water  column  float  array  a  large  distance  out  of  the  triangle  formed 
by  the  bottomed  floats  does  not  appreciably  affect  either  filters’  performance. 

•  When  the  mid-water  column  array  remains  within  the  triangle  formed  by  the  bottomed 
floats,  losing  one  of  the  bottomed  floats  does  not  appreciably  affect  either  filters’  per¬ 
formance. 

•  Deploying  the  mid-water  column  floats  in  a  vertical  line  array  increases  the  error  in 
the  Kalman  and  least  square  filter  estimates  by  30  and  50%,  respectively.  When  the 
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floats  are  stacked  vertically,  tlieir  interfloat  ranges  contain  little  information  about 
horizontal  displacements,  an<l  the  filters  are  forced  to  rely  exclusively  on  ranges 
between  the  mid-water  column  and  bottomed  floats  to  determine  the  x  and  y  coordi¬ 
nates  of  the  drifting  floats.  'I'he  Kalman  filter  fairs  somewhat  better  than  the  least 
squares  filter  because  of  its  ability  to  track  the  floats’  motion. 


6.  Summary 

We  have  evaluated  a  series  of  programs  which  proce.ss  Swallow  float  travel  time  meas¬ 
urements  to  produce  float  position  estimates.  Simulations  were  run  to  compare  the  perfor¬ 
mances  of  a  Kalman  filter  and  a  least  squares  estimator  under  several  deployment  geometries 
and  noise  conditions.  The  Kalman  filter  performed  better  than  the  least  squares  filter  under 
every  simulation  in  which  the  floats  experienced  small  random  accelerations  (low  process  noise). 
When  the  floats  experience  greater  accelerations  between  measurements,  the  Kalman  filter  esti¬ 
mate  error  increases  and  approaches  that  of  the  least  squares  estimator. 

The  Kalman  filter  functions  better  than  the  least  squares  filter  in  the  presence  of  high 
measurement  noise  when  the  process  noise  is  low.  This  is  because  of  its  ability  to  track  float 
motion  and  effectively  smooth  noisy  measurements.  When  both  the  process  and  measurement 
noise  are  increased,  the  Kalman  filter  estimation  error  approaches  that  of  the  least  squares 
filter. 

Both  filters  were  found  to  be  generally  insensitive  to  changes  in  deployment  geometry 
which  might  have  been  expected  to  increa.so  estimation  error.  This  is  in  part  because  our  simu¬ 
lation  did  not  include  error  in  the  average  sound  speed  estimate,  which  would  magnify  errors  in 
long  travel  times.  One  geometry,  in  which  the  mid-water  column  floats  were  reconfigured  from 
a  planar  array  to  a  vertical  line  array,  did  increase  the  estimation  error  in  both  filters, 
apparently  because  it  reduces  information  about  horizontal  float  displacements. 
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Float  4,  86  deployment;  range  from  floot  5 


Figure  5.  Simulating  Swallow  float  travel  time  measurements. 


Figure  6.  Swallow  float  localization  simulations. 
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Figure  7.  Effect  of  Process  and  Measurement  Noise 
on  Filter  Performance 
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Figure  8.  Effect  of  Deployment  Geometry 
on  Filter  Performance 
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Three-Dimensional  Beamforming  Results 


Power  Density 
99102 


Beamformer  Output  vs.  Look  Direction  fingles 
Record  1802,  9  Hz. 


129 


Elevation  Angle 


•I* 


ler  Output  vs.  Look  Direction  Rngles 
.802,  10.3  Hz. 


Elevation  Angle 


Power  Density 

90  93 


Beamformer  Output  vs.  LooU  Direction  Angles 
Record  1802,  12.2  Hz. 


Rzimuth  Rngle,  Degrees  CCW  from  X  Axis 


131 


Elevotxon  Rngle 


Power  Density,  dB 


Beomformer  Output  vs.  Look  Direction  Rngl 
Record  1844,  9  Hz. 


160  170  160  1^0  2^  2l0  : 

Azimuth  Angle,  Degrees  CCW  from  X  Axis 


Elevation  Angle 


Power 


Beamformer  Ouipui  vs.  Look  Direction  Angles 
Record  1844,  10.3  Hz. 


Beornfarmer  Ouipui  vs.  Look  Direction  Rngl 
Record  1844,  12.2  Hz. 


Azimuth  Angle,  Degrees  CCW  from  X  Axis 


Power  Densiiy,  dB 
102 


Beamformer  Oulpul  vs.  Look  Direction  Rngles 
Record  1844,  13.5  Hz. 


Azimuth  Angle,  Degrees  CCW  from  X  Axis 


135 


Elevation  Angle 


APPENDIX  E 


Geophone  Cross-Correlation 


IN 


200-175-150-125-100  -75  -50  -25  0  25  50  75  100  125  150  175  200 


ss- 


Cross-Correlation  of  9,  10.3,  12.2,  13.5,  16,  S(  22.5  Hz,  Record  1820 
Chomel  Impulse  Response,  Floats  9a,  24  Bin  Mask 


25  50  75  100  125  150  175  200 


Chiaf  of  NkTsl  Rmarch 
0«p»rtawat  of  iho  Notj 
Arliacloa,  Virciaik  3817-UaO 
Codo  1%  133(2),  lOaC 
in.  ItZ  US, 

U2SPO,  43S-AC,  460 

ONRDET 
NSTL  SlaiioB 

Bor  Si.  UmIo,  MiaiMippi  3S629-U04 
Codo  nZ  1121,  1122CS,  422CB, 
1122PO,  112600 

Diroctor 

OBco  of  Nonl  Ronwch 
Broach  OBco 
1090  Loot  Orooo  Strwt 
Poudoao,  Coliforaio  91106-2486 

Coconwador 

Novol  Soo  Sroioma  Coounoad 
Woohiacloo,  D.  C.  20363 
Cod*  63,  fSR,  63R-23 

Commoadlag  OBcor 
NotoI  Ocoaa  RoMorch  oad 
DoTolopoMat  Activilr 

NSTL  Siotioa 

Bay.  SI.  Looia,  Mlaiiiippi  3S638-6004 
Cod*  100,  lia  m  330, 

340,360,  360,  600 

Cc«uaaadw 

U.S.  Novol  Oeoaaoarophie  OBco 
NSTL  Slolioa 

Bor  Si.  Looio,  Mlanaiippi  38622-6004 
Alloa:  BUI  Jofect 

Aaiiloal  Socrolorr  of  tho  Navy 
(Roaoorch  Eaclaooriaf  A  Sytioiai) 
DcponsMBl  of  iho  Navy 
Woahiacloa,  0.  C.  20360 

D«fta«o  Advoacod  Roi.  Proj.  Afiacy 
TTO  -  Toclicol  Tochaolopy  OBco 
1400  Wiboa  Bovlovord 
ArUagloo,  Vlrdaia  2230^3308 
Alloa:  Copl.  Kirk  Evooi 

Nolioaol  Ocoaaic  A  Almoophoric 
Admiaiilraiioa 
Occoa  EasiBoorias  OBco 
6001  Ezocaliv*  Boalovard 
Rockvillo,  Marylaad  20863 

Cammoador 

Spoco  oad  Novol  Woiforo 
Syiiom*  Commaad 
Waohiafloa,  D.  C.  20360 
Codo  PDW.124,  32QA 

Commoador 

Novol  Ship  Roo.  A  Dov.  Coalor 
Bcihcfdo.  Marylaad  20084 

Exocuiivo  Socroiory 
Novol  Slodio*  Board 
Nolioaol  Acodomy  of  Scioaca 
2101  Cdoitilolioa  Avoauo,  N.W. 
Woohiailoa,  D.C.  20418 

Diroclor 

Siraiofic  SyoloiiM  Proj.  Ofc.  (PM-1) 
Doponiomi  of  Iho  Navy 
Woohiaploa,  0.  C.  30381 
Codo  NSP-2023 


ONR/MPL  GENERAL  DISTRmUTION  LIST 


Commoador 

Novol  Occoa  Sycicmi  Coalor 
Soo  Diofo.  Coliforoto  P21S2 
Codo  00,  01.  16.  94.  631 
6301.  Tl.  72 

Commoador 

Submoriao  Dovclopnoai  Group  ONE 

Flooi  Pool  OBco 

Soa  Dioco,  Coliforaio  92163 

Commoadias  OBcor 
CivU  EailBooriag  Lohoroiory 
Naval  Caoolractioo  Bollolio*  Coalor 
Port  Huoaomo,  Coliforaio  93043 
Codo  L40,  L42 

Commoadlag  OBcor 
Novol  Uadonrotor  Synomi  Coalor 
Noorport.  Rhodo  tolaad  20844 
Alloa:  Dr.  K.  A.  Komp,  Codo  8213 

OBcor  ia  Charg* 

Novol  Uodonvoior  Syrtomo  Coaler 
Now  Loadoa  Lohoroiory 
Now  Loadoa.  Coaaociicai 
Codo  900.  906.  910,  930.  960 

Diroctor  of  Roooorch 
U.S.  Naval  Roaaarch  Lohoroiory 
Woohiagloa,  0.  C.  20376 
Codo  2830.  2827.  6000,  5100,  6800 

Commoador 

Naval  Sarfoc*  Combol  Syotomo  Coaler 
While  Oak 

SUvor  Spriag.  Moiyload  20910 

Commoadlag  OBcor 

Novol  Coaoial  Syotomo  Loboroiory 

Poaoma  City,  Florida  32401 

STOtAC 

BailoUo  Columbuo  Lohorolorieo 
306  Kiag  Avoaoo 
Colambo*.  Ohio  43201 

Commoador 

Novol  Air  Syrtom*  Commaad 
Woohiagloa,  0.  C.  20381 
Codo  370 

Commaadiag  OBcor 
U.S.  Naval  Air  Dowlogmoai  Coalor 
Alioatioa:  Jim  Howard 
Warmiaiator,  PoaaiytTaaio  18974 

Diroclor 

Defoaoo  Docamoatotiaa  Coalor 
(TIMA),  Comoroa  Sloiioa 
6010  Dako  Street 
Alaxaadria,  Virgiaia  22314 

lanilolo  for  Dofoom  Aaalyrea 
1801  North  Boaagard  Siroot 
Arliagtoa,  Virgiaia  2231 1 

Saporiaioadoat 

U  S.  Naval  Pootpodaoio  School 
Moaloroy.  Califoraia  93940 

Chief  Scioaliot 

Navy  Uadorwaior  Soaad  Roforeaco  Div. 
US.  Naval  Racoarch  Laborolory 
P.O.  Boa  8337 
Orloado,  Florida  32808 


Supremo  Allied  Commoador 
U.S.  Alloalic  Fleoi 
ASW  Roooorch  Coaler.  APO 
Now  York.  Now  York  09016 
Via:  ONR  100  M.  CNO  OP092D1. 
Socrolarial  of  Miliiory, 
lafonnalioa  CoalroL  Commilteo 

Diroctor 

In-titato  of  Mariao  Sciooce 
Uaivomiy  of luka 
Foirbaak*,  Alaoko  99701 

Diroclor 

Applied  Phyoico  Lohoroiory 
Johao  Hopkia*  Uoivorriiy 
Joha*  Hopkia*  Rood 
LouroL  Marylaad  20810 
Alloa:  J.  R.  Aotlia 

Diroctor 

CoUog*  of  EagiBeeriag 
Doporimoai  of  Occoa  Eagiaoeriag 
Florida  Atlaatic  Uaiversiiy 
Boca  Raloo,  Florida  33431 

Diroctor 

Maria*  Roaaarch  Laboroiorie* 
c/o  hdaria*  Stadia*  Coalor 
Uanaialy  of  WiaeoaaiB 
Modboo,  Wlacaaria  63706 

Dlractor 

Applied  Roaaarch  Laboratory 
Poaaoylvaoia  Slot*  Uaivenrity 
P.O.  Boa  30 

Slot*  CoUogi,  Peaaqrivoaio  16802 
Diroctor 

Applied  Phyric*  Loboroiory 
Uaivoraity  of  Wubiagioa 
1013  Eaol  40lh  Siroot 
Soaltl*.  Woohiagloa  98106 

Diroctor 

Tho  Uaiv.  af  Toaa*  at  Auoiia 
Applied  Roaaarch  Loborotoiy 
P.O.  Baa  8029 
Auotia.  Toaa*  78712 

Diroctor 

Lomoat-Dohorty  Goologicol  Obrorvotory 
Torroy  CUI 

Polioadoa,  Now  York  10964 
Diroctor 

Wood*  Hole  Oceaaographic  laotitutioD 
Wood*  Helo,  Maaiacho**ita  02643 

Diroctor 

Ian.  of  Ocoaa  Scioaco  Eagiaooriag 
Colhelic  Uaivorrity  of  America 
Waakiagtoa,  D.C.  20017 

Nolioaol  Scioac*  Foaodolioa 
Waakiagtoa.  D.  C.  20660 

OBc*  of  Naval  Reooorch 
Rooidoat  Roproaoalotiv* 
c/o  Uaiv.  of  Coliforaio,  Soo  Diego 
Nkil  Cod*  Q033 
La  Jolla,  Califoraia  92093 

Uaiveraity  of  CoUforaio,  Soa  Diego 
Maria*  Phyoicol  Laboratory 
Broach  OBco 
La  Jolla,  Califoraia  92093 


May  1986 


